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I.  INTRODUCTION  AND  SUMMARY 


1.1  INTRODUCTION 

During  the  past  six  months  research  efforts  at  S-CUBED  under 
contract  F08606-80-C-0016  have  focused  on  automatic  discrimination 
and  surface  wave  path  corrections.  This  report  summarizes  the  major 
accomplishments  achieved  to  date  on  these  two  topics. 

The  objective  of  the  automatic  discrimination  effort  is  to 
design,  code  and  test  a  computer  program  which  analyzes  seismograms 
automatically  and  then  classifies  them  as  either  explosion-like  or 
earthquake-like,  with  associated  misclassification  probabilities. 
The  coding  and  implementation  of  the  program  for  automatic 
discrimination  is  following  the  approach  used  previously  for 
"Automatic  Seismic  Signal  Processing  Research"  (VT/0704).  The 
on-line  code  consists  of  two  modules.  One  performs  feature 
selection  and  produces  as  its  output  a  file  of  discriminants  (or 
features)  for  each  input  seismogram.  The  second  module 
statistically  combines  the  features  at  each  station  to  perform 
individual  station  classifications  and  then  polls  all  stations  to 
form  a  network  consensus  as  to  the  character  of  an  event. 

In  order  to  test  the  Automatic  Discrimination  code,  we  have 
been  assembling  an  expanded  collection  of  Area  of  Interest  (AI) 
seismograms  and  processing  them  in  a  manner  similar  to  the  original 
VELA  Seismological  Center  (VSC)  sponsored  AI  experiment.  The 
objective  of  this  evaluation  is  to  determine  the  effectiveness  of 
the  various  discriminants  with  respect  to  availability  and  type  of 
data,  source  size,  propagation  path  and  station  properties. 

The  objective  of  the  surface  wave  path  corrections  effort  is 
to  develop  computer  software  for  the  VSC  seismic  system  to  invert 
observed  surface  waves  for  path  structure  and  attenuation,  and  to 
provide  Green's  functions  for  moment  tensor  inversion  for  selected 
paths.  This  task  has  been  largely  completed  and  an  S-CUBED  Topical 
Report,  written  by  Stevens,  et.  aU  (1982)  and  submitted  to  VSC 
during  this  reporting  period,  gives  a  fairly  complete  description  of 


the  surface  wave  analysis  package.  The  report  also  serves  as  a 
user's  guide  for  the  software  package  which  was  made  operational  at. 
the  Seismic  Research  Center  (SRC)  during  the  past  six  months.  The 
reader  is  referred  to  appendices  to  the  report  by  Stevens,  et  _al_. 
(1982)  which  contains  file  formats,  sign  conventions,  an  explanation 
of  some  of  the  algorithms  used  in  the  programs,  and  definitions  of 
output  quantities. 

In  brief,  the  method  adopted  in  the  surface  wave  analysis 
package  uses  dispersion  to  obtain  the  travel  path  structure  which, 
in  turn,  is  used  to  estimate  the  moment  of  an  explosion  source  and 
path  attenuation.  The  technique  was  originally  used  by  Bache,  et 
al.  (1978)  to  estimate  the  moment  and  yield  of  NTS  events.  The 
basic  algorithms  used  in  that  study  were  subsequently  improved, 
extended,  and  simplified  so  that  path  corrections  may  be  obtained 
much  more  easily  in  an  interactive  environment. 

1.2  SUMMARY 

The  major  result  of  the  surface  wave  studies  described  in  this 
report  is  the  development  of  an  interactive  computer  program  for 
processing  long-period  seismic  data  and  its  application  to  the 
analysis  of  SRO  recordings  of  East  Kazakh  events.  The  analysis 
interprets  an  observed  vertical  component  long-period  seismogram  in 
terms  of  the  following  quantities. 

1.  The  dispersion  of  the  group  and  phase  velocities. 

2.  The  plane  layered  structure  (b,  Q)  which  best  explains 
the  dispersion. 

3.  The  moment  of  the  source. 

4.  The  matched  filter  which  best  removes  the  dispersion. 

In  addition,  a  synthetic  seismogram  module  permits  the  calculation 
and  display  of  the  seismogram  caused  by  an  explosion  source  in  a 
plane  layered  earth  model. 

The  conclusion  from  the  SRO  surface  wave  study  is  that,  with 
respect  to  East  Kazakh  KONO,  SH 10  are  excellent  stations,  MAJO,  GRFO 


are  good  stations,  ANTO,  CHTO  are  noisy  with  more  multipathing  than 
the  first  four,  and  KAAO  (distance  18  degrees)  is  very  difficult  to 
process.  CHTO  seismograms  have  a  peculiar  spectrum,  with  a  large 
dip  at  about  20  seconds  period. 

The  major  result  of  the  automatic  discrimination  studies 
described  in  this  report  is  the  development  of  a  new  multivariate 
body  wave  discriminant  based  upon  numerical  source  simulations. 
Source  identification  using  this  method  is  called  "deterministic" 
discrimination,  because  the  definition  of  the  discriminant  depends 
entirely  upon  deterministic  body  wave  modeling.  Application  of  the 
discriminant  to  a  modest  data  set  shows  good  results.  Furthermore, 
the  method  of  deterministic  discrimination  can  be  used  to  estimate 
the  body  wave  attenuation  parameter,  t*,  and  we  present  t*  result 
for  SRO  recordings  of  NTS  and  East  Kazakh  sources. 

The  conclusion  from  the  automatic  discrimination  study  is  that 
deterministic  source  models  can  effectively  serve  as  training  data 
in  a  multivariate  discriminant.  Furthermore,  the  deterministic 
discriminant  can  be  viewed  as  a  way  of  estimating  the  excess  high 
frequency  attenuation  (t*)  that  occurs  when  waves  propagate  in  the 
earth  instead  of  the  computer.  Measurements  of  t*  using  the  new 
technique  agree  with  independent  studies  of  the  same  data. 
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II.  SURFACE  WAVE  ANALYSIS  PACKAGE 


The  following  is  a  description  of  the  S-CUBED  Surface  Wave 

Analysis  Package  which  is  designed  to  obtain  path  corrections  from 

obseved  seismograms.  There  are  three  main  programs  in  the  package. 

The  first  —  TELVEL,  uses  the  method  of  phase-matched  filters 

(Herrin  and  Goforth,  1977)  to  recover  surface  wave  phase  (C)  and 

group  (U)  velocities.  The  second  —  INVERT,  inverts  the  phase  and 

group  velocities  to  obtain  the  earth  structure  for  the  path,  and 

obtains  an  approximate  Q  (attenuation)  structure  and  moment  by 

comparing  synthetic  and  observed  spectral  amplitudes.  The  third  — 

SYNSRF,  uses  a  variety  of  methods  (Takeuchi  and  Saito,  1972;  Schwab 

and  Knopoff,  1970;  Harkrider,  1964,  1970)  modified  for  high 

frequency  stability,  to  synthesize  surface  waves. 

* 

All  programs  are  interactive  and  self-contained.  Every  effort 
has  been  made  to  make  them  as  easy  to  use  as  possible.  In  this 
report,  we  show  the  uses  of  the  programs  for  making  path  corrections 
by  going  through  an  example  in  detail  using  a  synthetic  seismogram, 
and  by  applying  the  method  to  Shagan  River  —  SRO  station  travel 
paths. 

2.1  WHAT  IS  A  PATH  CORRECTION? 

A  path  correction  is  a  Green's  function  for  a  given  source 
region  and  source-to-receiver  path.  Once  the  Green's  function  is 
known,  it  may  be  used  to  generate  synthetic  surface  waves  from  an 
arbitrary  source,  as  observed  at  the  receiver  point,  or  the  inverse 
problem  may  be  solved:  an  observed  seismogram  may  be  used  to 
estimate  the  strength  and  type  of  the  source,  usually  in  the  form  of 
a  moment  tensor. 

The  basic  procedure  used  to  make  a  path  correction  was 
outlined  in  an  earlier  report  (Wang,  et  aL,  1981).  Since  then  the 
procedure  has  been  improved  by  allowing  simultaneous  estimation  of  Q 
and  moment,  and  by  making  the  codes  interactive  and  compatible. 
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2.2  HOW  TO  MAKE  A  PATH  CORRECTION 


i 


I 


I 


To  make  a  path  correction,  simply  take  a  seismogram,  use 
program  TELVEL  to  obtain  the  Rayleigh  wave  phase  and  group 
velocities,  then  use  program  INVERT  to  find  the  average 
source-to-receiver  shear  velocity  (e)  and  Q  structure,  and  finally 
use  program  SYNSRF  to  generate  eigenfunctions  for  this  structure, 
and  to  compute  a  synthetic  seismogram  which  recembles  the  original 
seismogram.  The  scalar  moment  of  the  explosion  is  found  together 
with  the  estimate  of  Q,  and  the  eigenfunctions  can  be  used  to 
perform  a  moment  tensor  inversion  for  any  source  in  the  same  area  as 
the  original  explosion. 

In  practice,  finding  a  path  correction  is  more  than  a  routine 
operation,  and  there  are  a  number  of  ways  to  err  during  the  data 
processing.  The  purpose  of  thir  section  is  to  go  through  the 
procedure  in  some  detail  and  to  identify  the  traps  which  are  lurking 
to  confuse  the  user.  Fortunately,  most  of  the  errors  which  can  be 
made  will  show  up  somewhere  during  the  processing.  A  bad  set  of 
phase  and/or  group  velocities  may  not  produce  a  reasonable  earth 
model  when  inverted.  An  inaccurate  inversion  for  shear  velocity  may 
make  it  impossible  to  determine  a  reasonable  Q  structure.  If  the 
seismogram  obtained  at  the  end  of  the  procedure  is  a  good  match  to 
the  original  (except  for  noise  or  multipathing  effects,  etc.),  then 
chances  are  very  good  that  the  path  correction  has  been  accurately 
found. 

The  first  step  in  making  a  path  correction  is  to  collect  all 
of  the  information  pertaining  to  the  seismogram.  The  following 
items  are  required: 

1.  A  seismogram  sampled  at  evenly  spaced  points  in  REALIO 
format  demeaned,  detrended  and  tapered. 

2.  The  number  of  points  in  the  seismogram. 

3.  The  sampling  interval  for  the  seismogram. 

4.  The  time  delay  between  the  source  time  and  the  start  of 
the  seismogram. 
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5.  The  distance  from  the  source  to  the  receiver. 
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6.  The  instrument  response  expressed  as  a  ratio  of 
polynomials  also  in  REAL  10  format. 

It  is  very  important  to  have  an  accurate  estimate  of  source  to 
receiver  distance,  and  the  seismogram  start  time.  It  is  well  worth 
double  checking  these  numbers  since  everything  that  follows  will  be 
wrong  if  these  numbers  are  incorrect.  Of  course,  errors  of  a  few 
seconds  or  a  few  kilometers  may  riot  result  in  major  errors  in  the 
final  results,  but  it  is  the  abundant  typos  and  minute  errors  that 
can  be  disastrous.  The  best  way  to  be  sure  of  timing  is  to  compare 
several  seismograms  for  the  same  path  to  be  sure  the  times  are  all 
consistent. 

Several  small  utility  codes  are  available  which  allow  the 
conversion  of  any  seismograms  to  REAL  10  format  as  well  as 
detrending,  demeaning,  etc.,  and  which  put  instrument  responses  in 
the  proper  format  (see  appendices  in  topical  report  (Stevens,  et 
al.,  1982)  for  file  formats,  etc.).  Instrument  responses  may  be  in 
polynomial  form  or  may  be  tabulated  instrument  responses.  One 
additional  quantity  is  required  by  TELVEL  —  the  initial  phase  of 
the  source.  For  an  explosion  (vertical  component,  displacement, 
positive  up,  Rayuleigh  wave)  the  initial  phase  is  -  3»/4.  In  order 
to  illustrate  the  procedure  we  apply  the  entire  surface  wave 
analysis  package  to  a  synthetic  seismogram.  This  is  a  useful 
exercise  to  go  through  before  processing  data  for  a  new  area,  as  the 
results  help  to  identify  certain  problems  in  advance,  and  may  help 
to  Identify  the  correct  2*  branch  in  the  phase  velocity  analysis. 
If  an  approximate  structure  is  available  for  an  area,  it  should  be 
used  to  make  a  similar  test. 

A  synthetic  seismogram  was  constructed  (with  program  SYNSRF) 
using  the  structure  listed  and  shown  in  Figure  1.  The  synthetic  was 
computed  using  an  SRO  long-period  instrument  at  a  distance  of  3000 
kilometers  and  an  explosion  source  with  a  T*,  (the  long-period  limit 
of  the  reduced  velocity  potential)  of  one  cubic  meter.  The 
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RHO 
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Q 

noo.o 
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2500.0 
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3300.0 
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8200.0 
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3300.0 
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8000.0 
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3300.0 
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100000.0 
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4200.0 

3200.0 

90.0 
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8000.0 
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80000.0 

8500.0 

4700.0 

3500.0 

200.0 

Figure  1.  East  Kazakh  shear  velocity  structure  used  for  synthetic 
seismogram. 
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synthetic  seismogram  and  true  phase  and  group  velocity  dispersion 
curves  are  shown  in  Figure  2.  Figure  3  shows  a  form  summarizing  the 
inputs  and  outputs  to  the  programs.  It  is  useful  to  have  a  similar 
form  available  before  starting  the  procedure. 

2.2.1  Recovery  of  Phase  and  Group  Velocities  and  Spectral  Amplitudes 

With  this  information  in  hand,  the  next  step  is  to  run 

TELVEL.  The  first  phase  of  TELVEL  is  to  apply  a  set  of  narrow-band 

filters  to  the  seismogram  to  obtain  an  approximate  set  of  group 

velocities.  The  choice  of  filter  frequencies  and  filter  widths  j 

(specified  in  terms  of  the  filter  Q)  is  up  to  the  user.  From 

experience,  a  Q  of  10  provides  a  good  filter  width  for  most  data./ 

For  long-period  data,  such  as  SRO  data,  20  to  40  frequencies  from 

0.01  to  0.1  Hz  usually  produces  a  good  group  velocity  curve  over  the 
« 

reliable  frequency  range  of  data  (see  Figure  4).  Occasionally,  the 
program  will  encounter  difficulties  in  finding  a  good  group  velocity 
curve.  If  it  does,  try  going  back  and  using  more  or  fewer 
frequencies  or  a  lower  or  higher  Q.  Sometimes  the  program  will  find 
a  good  set  of  group  velocities  as  evidenced  by  the  peaks  on  the  j 

plot,  but  will  not  recognize  it  (no  curve  will  be  drawn  through 
them).  If  so,  they  can  be  entered  by  hand  using  the  edit  command. 

The  narrow  band  filter  is  an  important  step  in  finding  the  group  and 
phase  velocites.  Not  only  does  it  provide  an  initial  group  velocity 
estimate,  but  it  may  show  up  features  such  as  bifurcated  group 
velocity  curves  which  may  render  the  phase  matched  filter  unstable. 

The  purpose  of  TELVEL  is  to  improve  the  initial  estimate  of 
group  velocity  by  phase  matched  filtering  while  eliminating 
interfering  multipath  arrivals,  higher  modes  and  noise,  and  to  find 
the  phase  velocities.  The  phase  matched  filter  is  found  by 
integrating  the  group  delay  found  initially  by  narrow  band 
filtering.  When  the  filter  is  applied  to  the  original  seismogram, 
the  true  surface  wave  is  compressed  into  a  narrow  time  window.  The 
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VELOCITY  km/sec 


Source  Location:  East  Kazakh 

Receiver  Location:  Synthetic  a  3000  km 

Source  Time:  0  sec 

Receiver  Time:  600s ec 

Delay  Time:  600  sec 

Distance:  3000  km 

A 

Seismogram  File:  RSYN  *  EKSRO 
Number  of  Points:  901 
&T:  1.0  sec 

Instrument:  SRO-LP 
Initial  Phase  »  -  0.75 

Narrow  Band  Filters  Used 
NF,  FMIN,  FMAX,  Q  -  30,  .01,  .1,  10 

TELVEL  Output  File(s):  TEL  *  EKSRO 

Eigenfunction  File  for  Source  Region  £F  *  EKAZ2 
Instrument  Gain  (conversion  to  meters)  1. 
Estimated  Moment  ■  6.4  x  10**  / Gain  »  6.4  x  10** 
Estimated  f,*  ■  .980  /Gain  »  .980 

INVERT  Output  File(s):  INV  *  EKSRO 

SYNSRF 

Frequencies  Computed:  100,  .002,  .2 
Output  Files:  EF  *  EKSRO 
Estimated  M$:  0.04 

Comments:  Test  on  Synthetic  Seismogram 
PATH  CORRECTION  INFORMATION 

Figure  3.  Form  summarizing  inputs  and  outputs  for  Surface  Wave 
Analysis  Package. 
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Figure  4.  Group  velocity  curve  obtained  by  narrow  band  filtering  the  synthetic  seismogram  in  Figure 
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resulting  phase  matched  filter  output  is  plotted  by  TELVEL  to  allow 
windowing  of  the  time  series  (see  Figure  5).  It  is  important  to 
choose  a  time  window  large  enough  to  contain  most  of  the  energy  in 
the  main  arrival,  but  small  enough  to  eliminate  most  of  the  noise 
and  multiple  arrivals  from  the  filter.  The  best  way  to  choose  the 
time  limits  is  to  look  at  the  phase  matched  filter  output  and  then 
back  up  to  apply  the  appropriate  time  window.  For  long  period  data 
a  time  window  of  50-100  seconds  on  each  side  of  the  main  pulse  is 
usually  appropriate. 

After  the  phase  matched  filter  is  applied  to  the  seismogram, 
the  filter  phase  is  unwrapped  and  used  to  estimate  the  phase 
velocities  of  the  Rayleigh  wave.  Since  any  multiple  of  2ir  may  be 
added  to  the  phase,  there  is  an  uncertainty  introduced  into  the 
phase  velocity.  If  6  is  the  phase  of  the  surface  wave  (d  is  always 
<  0)  then  the  phase  velocity  is  given  by 

r  -  2wfr 

u  +'"2'irn 

where  r  is  the  distance,  tQ  is  the  initial  phase,  f  is  the 
frequency  and  n  is  an  undetermined  integer.  The  choice  of  the 
correct  value  of  n  is  often  difficult.  In  order  to  help  choose  the 
correct  value  of  n,  TELVEL  plots  the  phase  velocity  for  a  user 
specified  value  of  n  together  with  values  of  n  *  1  (see  Figure  6). 
The  following  facts  help  in  the  choice  of  n. 

1.  The  true  phase  velocity  approaches  a  constant  in  the  low 
frequency  limit.  If  the  phase  velocity  is  well 
determined,  values  of  n  which  are  too  small  cause  the 
phase  velocity  to  decrease  rapidly  at  low  frequencies, 
while  values  of  n  which  are  too  large  cause  the  phase 
velocity  to  increase  rapidly  at  low  frequencies.  In 
practice,  the  correct  phase  velocity  is  usually  the  first 
curve  to  turn  upwards  at  low  frequencies,  although  it  is 
sometimes  the  first  to  turn  downwards,  and  there  is 
sometimes  no  clear  choice. 
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ase  matched  filter  output  for  synthetic  seismogram.  Energy  is  compressed  to  within 
100  seconds  of  primary  arrival  time  based  on  initial  group  velocity  curve. 


OPTIONtM  FOR  HELP)  >L 
1M  FREO  PHftSE-W 


2.  The  true  phase  velocity  is  not  very  different  from  the 
group  velocity  at  low  frequencies,  and  is  always  greater 
than  the  group  velocity  at  the  lowest  accurately 
determined  frequency. 

3.  The  phase  velocity  is  approximately  equal  to  4  km/sec  in 
most  regions  of  the  world  at  frequencies  between  0.01  and 
0.02  Hz. 

If  the  choice  of  number  of  2¥*s  to  add  to  the  phase  is  not 
clear,  save  output  files  for  more  than  one  value  and  see  which  one 
leads  to  a  reasonable  inversion  model.  Also,  processing  of  another 
seismogram  for  the  same  path  may  lead  to  a  clearer  choice  of  phase 
velocity. 

An  iterative  procedure  is  used  to  refine  the  group  velocity 
estimate.  Usually  two  or  three  iterations  will  produce  stable  group 
velocity  estimates  which  do  not  change  from  one  iteration  to  the 
next  (see  Figure  7).  It  is  a  good  idea  to  save  the  phase  velocity 
plot,  and  group  velocity  plot  on  the  final  iteration  for  future 
reference.  The  choice  of  2*'s  to  be  added  to  the  phase  may  be 
easier  after  the  group  velocity  has  converged. 

The  final  decision  to  be  made  is  the  choice  of  freqnencies  to 
output.  These  do  not  need  to  be  the  same  as  the  initial  narrow  band 
filter  frequencies.  The  range  of  frequencies  should  be  large  enough 
to  Include  all  data  points  that  are  reliable,  but  small  enough  to 
exclude  unreliable  data  points.  Two  figures  obtained  in  the 
processing  help  to  identify  the  appropriate  range.  Points  on  the 
initial  narrow  band  filtered  group  velocity  plot  with  amplitudes 
less  than  10  to  20  per  cent  (indicated  by  1  and  2  on  the  plot)  of 
the  maximum  amplitude  are  not  likely  to  be  reliable.  On  the  final 
group  velocity  plot,  portions  of  the  curve  at  the  ends,  which  have 
not  stabilized,  are  not  reliable.  Portions  of  the  group  velocity 
curve  which  appear  odd  such  as  those  which  have  a  rapid  decrease  in 
velocity  or  an  unusually  high  velocity  at  low  frequency  or  rapid 
oscillations  at  high  frequency  should  not  be  used.  Periods  longer 
than  the  time  window  used  are  also  unreliable.  It  is  important  to 
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obtain  information  at  low  frequencies  since  these  frequencies 
determine  the  deep  structure  of  the  earth.  For  most  SRO  station  at 
distances  of  2000  through  5000  kilometers,  phase  and  group 
velocities  may  be  reliably  obtained  from  about  0.015  to  0.08  Hz.  A 
convenient  frequency  spacing  is  0.005  Hz.  Phase  and  group 
velocities  can  almost  always  be  obtained  down  to  a  minimum  frequency 
of  0.02  Hz. 

The  output  file  from  TELVEL  contains  the  estimated  phase  and 
group  velocities  and  the  instrument  corrected  spectral  amplitudes  of 
the  seismogram.  The  corrected  amplitudes  should  be  an  improvement 
over  the  Fourier  amplitudes,  since  they  are  obtained  from  the 
compressed  seismogram  after  the  removal  of  noise  and  interfering 
phases. 

The  TELVEL  output  file  also  contains  a  space  for  the  standard 
deviation  of  the  phase  and  group  velocity.  TELVEL  does  not  estimate 
errors  in  the  velocities,  so  it  simply  outputs  default  values  o> 
zero.  If  more  than  one  seismogram  for  the  same  path  are  processed, 
they  can  be  averaged  to  get  an  improved  estimate  fit  fhase  and  ^roup 
velocities  and  to  estimate  uncertainties  in  the  data.  A  small 
utility  program  (AVEDAT)  is  available  which  constructs  a  combined 
file  in  the  format  of  a  TELVEL  file. 

2.2.2  Inversion  for  Earth  Structure 

The  TELVEL  output  file  is  in  the  proper  format  to  be  used  by 
the  second  major  code  —  INVERT,  which  inverts  the  phase  and  group 
velocities  for  the  average  shear  velocity  structure  along  the  path, 
and  which  uses  the  amplitudes  to  estimate  the  Q  structure  of  the 
path. 

Program  INVERT  has  been  designed  to  require  a  minimum  of 
operator  input.  A  starting  model  is  automatically  generated  from 
the  input  data.  The  final  model  is  independent  of  the  starting 
model,  but  it  is  important  to  choose  a  suitable  range  of  depths 
appropriate  for  the  resolution  capability  of  the  data. 

Only  the  TELVEL  output  file  is  required  for  shear  velocity 
inversion.  One  other  file  is  needed  for  Q  inversion.  Assuming  that 

17 


the  structure  of  the  source  region  is  approximately  known 
independently  from  previous  investigations,  an  eigenfunction  file 
(generated  for  the  source  region  structure  by  program  SYNSRF)  will 
be  input  at  this  stage. 

INVERT  is  divided  into  three  main  sections  —  model 
generation,  inversion,  and  analysis.  Most  of  the  time  the  only 
command  necessary  in  the  model  generation  section  is  RUN  which 
causes  model  generation  to  performed.  In  the  model  generation 
section,  the  user  can  set  a  value  of  Poisson's  ratio  and  the 

coefficients  of  Birch's  law  if  desired.  The  inversion  code  only 
inverts  for  shear  velocity  while  compressional  velocity  and  density 
are  constrained  by  these  constants.  The  default  values  are 
appropriate  for  most  of  the  earth  and  rarely  need  to  be  changed. 
Two  features  which  are  sometimes  needed  are  the  ability  to  set  a 
discontinuity  at  a  particular  depth,  and  the  ability  to  set  the 

number  of  layers  used  in  the  inversion.  The  default  number  of 
layers  is  20.  A  smaller  number  of  layers  will  result  in  a  faster, 
but  less  accurate  inversion. 

In  many  cases  a  discontinuity  may  be  required  to  match  the 

data,  but  its  depth  is  not  usually  known  in  advance.  The  best  plan 
is  to  run  the  inversion  section  and  insert  a  discontinuity  where  the 
resulting  model  shows  a  strong  velocity  gradient. 

The  first  step  in  the  inversion  is  to  constrain  the  smoothness 
of  the  model  defined  by  the  number  of  degrees  of  freedom  (command 
OF)  allowed.  A  small  value  of  OF  produces  a  smooth  model  while  a 

higher  value  produces  a  better  data  fit.  A  good  value  to  start  with 
is  a  OF  of  4. 

The  command  INVERT  causes  the  first  iteration  of  the  inversion 
to  be  performed.  This  command  causes  partial  derivatives  ac/as  and 
au/ae  to  be  calculated  from  the  initial  (or  current)  model.  It 
causes  a  set  of  matrices  to  be  assembled  from  the  partial 
derivatives,  and  performs  a  full  singular  value  decomposition  of  the 
matrices.  Finally  it  computes  the  exact  nonlinear  phase  and  group 
velocities  for  the  estimated  model  with  the  OF  specified  initially. 
All  of  this  requires  quite  a  lot  of  computation  and  will  take  some 
computer  time. 


After  the  inversion  has  been  performed,  the  command  PLOT  will 
plot  the  shear  velocity  structure  while  the  command  PLOT  DATA  will 
plot  the  linearized  data  fit. 

The  command  OF  has  a  new  meaning  after  an  inversion  has  been 
performed.  A  major  advantage  of  the  matrix  decomposition  is  that 
new  models  for  different  values  of  DF  may  be  generated  without 
recomputing  the  matrices.  Thus  if  the  command  DF  4  was  given  before 
inversion,  the  command  DF  5  after  inversion  instantly  produces  the 
(linearized)  model  for  a  DF  of  5.  The  commands  PLOT  and  PLOT  DATA 
again  produce  the  model  and  data  fit  for  the  new  DF  (see  Figure  8). 
This  makes  it  very  easy  to  pick  the  best  value  of  DF.  A  second 
iteration  can  then  be  performed  using  this  model  as  the  new  starting 
model.  Usually  three  iterations  with  the  same  DF  are  sufficient  for 
convergence. 

An  obvious  question  is  how  to  choose  the  best  model.  The 

answer  is  to  use  a  DF  high  enough  to  get  a  good  data  fit,  but  small 
enough  to  prevent  oscillations  of  the  shear  velocity  model,  and  to 
yield  a  smooth  fit  to  tne  data.  Usually  a  DF  of  about  5  gives  a 

good  model.  With  good  quality  data  such  as  the  SRQ  data  a  DF  of  6 

may  be  used,  while  less  accurate  data  may  only  allow  a  DF  of  4  (DF 
need  not  be  an  integer). 

One  other  choice  must  often  be  made  during  inversion.  It  may 
be  necessary  to  include  a  discontinuity  in  the  model.  The  most 

common  place  a  discontinuity  is  needed  is  at  the  crust-mantle 

boundary.  The  proper  location  for  the  discontinuity  may  be 

determined  by  the  midpoint  of  a  sharp  gradient  which  appears  in  the 
shear  velocity  model.  Within  the  inversion  section  a  discontinuity 
may  be  added  by  specifying  the  number  of  the  layer  above  the  desired 
discontinuity,  or  it  may  be  specified  by  depth  by  returning  to  the 
model  generation  section  (this  is  unnecessary  unless  the  user 

requests  a  specific  depth  which  does  not  coincide  with  a  layer 

boundary).  If  higher  frequency  data  are  available,  as  it  may  be  for 
WWSSN  instruments  at  distances  less  than  1000  kilometers,  it  may  be 
necessary  to  include  a  shallow  discontinuity  in  the  upper  few 


FREQUENCY 


Figure  8.  Shear  velocity  structure  obtained  by  inverting  dispersion 
curves  produced  by  TELVEL  from  the  synthetic  seismogram, 
together  with  linearized  data  fit.  The  strong  gradient  from 
20  to  60  km  Indicates  a  possible  discontinuity.  This  model 
Is  shown  for  the  second  iteration  using  a  DF  of  5.0. 

(X)  Indicates  synthetic  data, (U)  indicates  estimated  group 
velocity,  (C)  indicates  estimated  phase  velocity. 
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kilometers.  It  may  be  difficult  to  define  the  proper  depth,  since 
the  phase  and  group  velocities  may  lose  accuracy  at  the  high 
frequencies  needed  to  resolve  shallow  depths.  When  the  presence  of 
this  layer  is  indicated,  it  is  often  essential.  The  NTS-Tucson 
seismograms,  for  example,  absolutely  cannot  be  matched  without  a 
shallow  low  velocity,  low  Q  layer  in  the  upper  two  kilometers 
(Bache,  T.  C.,  W.  L.  Rodi,  and  0.  G.  Harkrider  (1978)). 

After  the  last  iteration  has  been  performed,  give  the  command 
END  to  finalize  the  inversion  and  go  to  the  analysis  section.  Here 
you  can  make  final  plots  of  the  model,  and  the  true  nonlinear  data 
fit  (Figure  9).  You  can  also  print  the  data  fit,  variances  and 
spreads  for  the  final  model. 

One  more  calculation  is  required  to  obtain  the  complete  earth 
model  —  the  Q  structure  needs  to  be  found.  The  program  does  this 
by  computing  a  synthetic  seismogram  (actually  an  unattenuated 
synthetic  spectrum)  at  the  observed  frequencies  and  taking  the  ratio 
with  the  spectrum  of  the  filtered  data.  The  logarithm  of  the  ratio 
is  equal  to  a  constant  (the  moment)  plus  a  frequency  dependent 
attenuation  coefficient.  In  effect  the  level  of  the  spectral  ratio 
gives  the  moment,  while  the  shape  of  the  spectral  ratio  as  a 
function  of  frequency  gives  the  attenuation  coefficients  which  can 
be  inverted  for  the  Q  structure.  As  mentioned  before,  if  an 
eigenfunction  file  is  available  for  the  source  region,  the  program 
will  generate  the  synthetic  spectrum  using  these  eigenfunctions  for 
the  source  region  and  using  the  inverted  structure  for  the  path. 
The  program  uses  the  spectral  ratio  and  inverts  for  log  moment  and 
B/Q  simultaneously. 

The  spectral  amplitudes  are  never  known  as  accurately  as  the 
phase  and  group  velocities,  so  a  lower  DF  must  be  used  for  Q 
inversion.  The  minimum  possible  DF  is  2.  A  DF  of  2  will  produce  Q 
proportional  to  s  in  all  layers.  An  important  function  of  the  Q 
inversion  is  in  fact  to  smooth  the  spectral  ratios  and  the 
attenuation  coefficients.  It  is  usually  necessary  to  use  a  DF 
between  2  and  3.  Higher  DF's  result  in  unrealistic  Q  structures. 
Discontinuities  are  not  allowed  for  Q  inversion,  and  are  not  used 
even  if  they  were  used  in  the  velocity  inversion. 
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Figure  9 

► 


Final  model  obtained  using  a  DF  of  6.0  and  a  single  discontinuity 
(compare  with  Figure  1)  together  with  data  fit.  The  intermediate 
layer  from  30  to  50  km  has  been  split  by  the  single  discontinuity. 
The  dispersion  curves  are  an  excellent  fit  to  the  data. 
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Q  inversion  is  linear,  and  takes  little  computer  time. 
Matrices  are  formed  and  decomposed  on  the  first  iteration,  so  the 

•  best  Q  model  may  be  quickly  found.  Again  the  command  PLOT  plots  the 
Q  model  while  the  command  PLOT  DATA  plots  the  spectral  ratios  and 
corresponding  data  fit.  The  spectral  ratios  should  look  like  an 
upward  sloping  line  (see  Figure  10).  If  the  data  points  do  not  have 

•  an  upward  slope,  it  means  the  attenuation  coefficients  decrease  with 
frequency  and  no  Q  structure  will  fit  the  data.  If  this  happens, 
there  is  probably  something  wrong  with  the  inversion,  or  with  the 
data. 

®  When  the  Q  inversion  is  performed,  the  estimated  moment  and  too 

of  the  explosion  are  printed  out.  Save  these  numbers,  since  they 
are  not  printed  out  again. 

^  Since  the  Q  inversion  produces  a  smoothed  a/Q  model,  it  will 

not  reveal  any  true  discontinuities  in  Q  in  the  earth.  In 

particular,  it  will  not  usually  produce  a  very  low  Q  at  the 
surface.  If  high  frequency  data  have  been  used,  it  may  be  necessary 
to  add  a  low  Q  layer  near  the  surface  after  the  inversion  to  reduce 
high  frequency  ringing. 

Once  the  inversion  is  completed  the  final  step  is  to  output  a 
structure  file.  The  structure  file  is  a  formatted  file,  which  is  in 

•  the  format  of  an  input  file  to  SYNSRF. 

2.2.3  Generation  of  Synthetic  Seismograms  and  Excitation  Functions 

The  third  stage  of  the  path  correction  procedure  is  to  use  the 

**  program  SYNSRF  to  generate  an  eigenfunction  file,  and  to  make  a 

synthetic  seismogram  as  the  final  test  of  the  procedure. 

SYNSRF  is  divided  into  four  main  sections.  The  first  section 

I  requires  input  of  a  model  (structure)  file  and  a  choice  of 

frequencies  for  the  calculation  of  phase  velocities.  The  second 
section  calculates  group  velocities,  eigenfunctions  and  related 
parameters.  The  third  section  calculates  a  synthetic  seismogram. 

^  The  fourth  section  allows  prints  and  plots  of  the  seismogram, 

spectra  and  dispersion  curves. 


fcv:-; 
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FREQUENCY 


Figure  10.  Q  structure  and  data  fit  obtained  using  a  DF  of  2.5.  The  Q 
model  is  a  smoothed  version  of  the  true  Q  model  (Figure  1) 
while  the  data  Is  the  ratio  of  synthetic  (unattenuated)  to 
observed  spectra.  The  moment  is  obtained  simultaneously.  F 
this  example  the  moment  differed  by  2%  from  the  true  moment. 


SYNSRF  is  designed  to  be  easy  to  use  and  to  make  it  easy  to 
correct  mistakes.  If  you  run  one  of  the  sections  and  use  an 
incorrect  number,  or  forget  to  output  a  file,  simply  back  up,  make 
the  cnange  and  run  this  section  again.  None  of  the  sections  is  very 
time  consuming  and  all  previous  inputs  are  remembered  (except  for 
output  files,  since  they  are  closed  after  each  section  is  run). 

The  essential  input  in  part  1  of  SYNSRF  is  the  model  file  (the 
output  file  from  INVERT).  The  default  set  of  frequencies  is  100 
frequencies  evenly  spaced  from  .002  to  .2  Hz.  In  most  cases,  this 
will  be  a  good  range  of  frequencies.  For  observations  at  close 
range  (  <  1000  km)  with  an  instrument  that  does  not  filter  out 
higher  frequencies,  it  may  be  necessary  to  use  frequencies  from  .005 
to  .5  Hz  instead.  The  program  is  currently  dimensioned  to  allow  100 
frequencies  and  three  modes.  After  the  model  file  is  read  in  in 
Part  1,  the  program  looks  through  the  shear  velocities  to  estimate 
root  search  parameters  (minimum  and  maximum  phase  velocities  and  a 
step  size).  The  default  values  are  almost  always  adequate, 
especially  for  the  fundamental  mode.  If  the  step  size  chosen  is  too 

large,  a  root  may  be  missed.  This  will  show  up  as  a  glitch  in  the 

dispersion  curve. 

The  command  RUN  causes  dispersion  curves  to  be  calculated  and 
leaves  the  program  in  Part  2.  The  only  inputs  needed  in  Part  2  are 
the  source  depth  and  a  file  name  for  the  eigenfunction  output  file. 
The  default  source  depth  is  1  kilometer.  If  the  structure  being 
used  is  the  path  structure,  then  the  source  depth  is  irrelevant.  It 
is  necessary  to  specify  an  eigenfunction  file  name  here.  This  file 
will  be  used  in  two  ways.  It  will  be  read  back  in  Part  3  if  a 

separate  source  region  is  used,  and  it  contains  all  of  the 

information  needed  to  construct  the  Green's  function  that  is  the 
final  path  correction. 

The  command  RUN  now  causes  the  eigenfunctions  to  be  generated 
and  leaves  the  program  in  Part  3.  You  can  make  a  plot  of  the 
dispersion  curves  (phase  and  group  velocities)  here  to  compare  with 
the  observations.  Make  this  plot  before  reading  in  a  new  source 


region  eigenfunction  file,  since  the  velocities  plotted  are  the 
source  region  velocities. 

Several  quantities  must  be  input  here  to  make  a  synthetic 
seismogram.  If  you  are  using  separate  source  and  path  regions,  an 
eigenfunction  file  for  the  source  region  and  the  eigenfunction  file 
for  the  path  must  be  read  in.  An  instrument  response  must  be  read 
in,  preferably  using  the  polynomial  coefficients  that  were  in 
TELVEL.  The  distance  and  time  window  must  be  specified  here.  After 
the  distance  is  set,  the  time  window  is  estimated  by  the  program 
using  the  group  velocities.  For  the  purpose  of  making  a  path 
correction,  however,  it  may  be  preferable  to  use  the  time  window  of 
the  original  seismogram  to  make  the  comparison  easier. 

The  command  RUN  now  causes  a  seismogram  for  an  explosion  with 
a  'I'ao  of  one  cubic  meter  to  be  generated,  and  leaves  the  program  in 
Part  4.  Then  plots  of  the  seismogram  and  spectrum  may  be  made,  a 
run  summary  may  be  obtained,  and  surface  wave  magnitudes  may  be 
computed . 

The  seismogram  plotted  here  should  look  like  the  original 
seismogram  (see  Figure  11)  after  removal  of  noise  and  multipathinq. 
The  amplitude  should  differ  from  the  original  amplitude  by  a  factor 
approximately  equal  to  the  ¥«,  estimated  in  the  inversion  calculation. 

The  last  step  in  the  procedure  is  to  transform  the 
eigenfunction  file  into  a  form  usable  by  a  moment  tensor  inversion 
program.  A  short  program  EXCITE  turns  the  eigenfunction  file  into  a 
formatted  file  using  the  excitation  function  notation  of  Kanamori 
and  Given  (1981).  A  more  detailed  explanation  of  this  approach  is 
given  in  Appendix  4  of  the  report  by  Stevens,  et_aj_.  (1982). 

2.3  SHAGAN  RIVER-SRO  PATH  CORRECTIONS 

In  this  section  we  describe,  in  brief,  the  processing  of 
explosion  seismograms  recorded  at  seven  SRO  stations.  For  each  path 
there  is  a  description  of  the  analysis  including  an  estimate  of 
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confidence  for  each  path  correction,  plus  five  figures.  The  first 
figure  is  a  narrow  band  filter  estimate  of  the  group  velocity  for 
one  seismogram  along  each  path.  The  NBF  estimate  often  shows  up 

problems  such  as  extreme  multipathing,  high  noise  level,  or 

bifurcated  group  velocity  curves  which  may  cause  difficulty  in  the 
data  analysis.  The  second  figure  is  a  plot  of  the  data  fit,  solid 

lines  being  the  calculated  phase  and  group  velocity  curves  from  our 

final  structure  while  letters  U  (group  velocity)  and  C  (phase 
velocity)  are  observed  data  points.  The  third  figure  is  a  plot  of 
shear  velocity  versus  depth  for  the  structure  obtained  by  inversion, 
the  fourth  is  a  plot  of  Q  versus  depth,  while  the  last  figure  is  a 
listing  of  velocities,  densities,  obtained  for  each  model.  The 
three  events  processed  are  presumed  explosions  at  the  Shagan  River 
test  site  which  occurred  on  2  December  1979  (No.  318),  29  November 
1978  (No.  312),  and  23  June  1979  (No.  313)  according  to  the  U.  S. 
Geological  Survey  Preliminary  Determination  of  Epicenters  Bulletin. 
Event  312  is  apparently  a  double  event,  in  which  an  mb  6.0 
followed  an  5.3  by  4.0  seconds.  We  used  the  parameters 

appropriate  for  the  larger  event  in  the  processing,  and  could  see  no 
sign  of  the  smaller  event.  The  seismograms  are  shown  in  Figures  12 
through  14. 

In  addition  to  obtaining  a  velocity  and  Q  structure  for  each 
path,  we  also  obtain  an  estimate  of  moment  for  each  path.  By 
assuming  aknown  congressional  velocity  and  density  for  the  source 
region  (we  used  a  -  5000  m/sec,  p  »  2100  kg/m^),  the  result  may  be 
expressed  in  terms  of  'f*  using  the  relation  Mx  »  4wp«  f*.  The 
results  are  listed  in  Table  1. 

The  last  station  is  considerably  less  reliable  than  the  first 
six  for  reasons  that  will  be  explained  with  the  data.  It  is 
Interesting  to  note  that  there  seems  to  be  a  correlation  between  the 
estimated  foo  and  the  quality  of  the  data.  This  may  be  an  indication 
of  energy  lost  due  to  scattering  and  non-plane  layered  effects  along 
the  travel  path.  It  does  not  seem  to  result  from  application  of  the 
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Figure  12.  Seismograms  recorded  at  SRO  stations  for  Shagan  River 

explosion  number  318,  December  2,  1979.  Seismograms  were 
processed  to  obtain  path  corrections  for  stations  KAAO, 

SHIO,  ANTO,  CHTO,  KONO,  GRFO  and  MAJO.  Before  processing, 

< $  a  reduced  time  window  was  chosen  to  eliminate  early  arrivals 

and  glitches  and  the  seismograms  were  demeaned  and  detrended. 
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Figure  13. 


Seismograms  recorded  at  SRO  stations  for 
explosion  number  312,  November  29,  1978. 
processed  to  obtain  path  corrections  for 
ANTO,  SHIO,  KONO,  and  GRFO. 


Shagan  River 
Seismograms  were 
stations  CHTO,  MAJO, 
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Seismograms  recorded  at  SRO  stations  for  Shagan  River 
explosion  number  313,  June  23,  1979.  Seismograms  were 
processed  to  obtain  path  corrections  for  stations  CHTO,  MAJO, 
SHIO,  and  KONO.  Processing  of  stations  KAAO  and  ANTO  was 
attempted,  but  excessive  noise  and  multipathing  made  it 
difficult  to  obtain  stable  group  velocity  curves  at  those 
stations. 


id 


q  Figure  14. 
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Table  1 
ESTIMATED  V*, 
(MOMENT  -r  6.60  x  1011) 


Station 

Oi stance 

Event  No.  318 

Event  No.  312 

Event  No.  313 

KONO 

4370 

15100 

9980 

15000 

SHIO 

2926 

13400 

9760 

13400 

MAJO 

4918 

13000 

5100 

8100 

GRFO 

4700 

11800 

6700 

— 

ANTO 

3739 

10300 

5440 

— 

CHTO 

3890 

7200 

5000 

6700 

KAAO 

1885  * 

5700 

— 

phase  matched  filter  time  window  since  synthetic  seismograms  made 
using  the  estimated  values  of  ¥«,  are  not  consistently  lower  in 
amplitude  than  the  observed  seismograms  based  upon  time  domain 
measurements  (see  Table  2).  Of  course,  the  estimates  of  ¥«>  also 
depend  on  the  accuracy  of  the  instrument  gain.  These  seismograms 
were  provided  by  Data  Services  at  the  Seismic  Research  Center  (SRC) 
after  conversion  to  nanometers  (at  20  seconds). 

In  general,  the  processing  resulted  in  reasonable  velocity 
structures  and  Q  structures  for  most  of  the  stations  and  should 
provide  a  good  estimate  of  the  path  structure  and  attenuation.  The 
results  are  shown  in  Figure  15  through  21. 


Table  2 

TIME  DOMAIN  AMPLITUDE  COMPARISON 


STATION  OBSERVED/SYNTHETIC  AMPLITUDES  (NM) 


Event  No.  318 

Event  No.  312 

Event  No.  313 

KONO 

765/725 

460/479 

736/720 

SHIO 

606/600 

414/437 

576/600 

MAJO 

371/378 

134/148 

248/236 

GRFO 

361/484 

240/275 

— 

ANTO 

192/156 

94/83 

143/ - 

CHTO 

174/249 

132/173 

178/232 

KAAO 

613/305 

_ 

282/ — 

Figures  15.1  through  15.5  are  on  the  following  pages. 


Figure  15.  Path  1:  SHAGAN-KQNO 

Distance:  4372  km 
Azimuth:  72.6* 

Instrument:  ASRO-LP 

Events  Processed:  318,  312,  313 

Description:  This  is  an  excellent  station.  A  very 

clean  group  velocity  curve  was  obtained  with  no  evidence 
of  multipathing  or  any  observational  problems. 
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SHAGAN-KONO  STRUCTURE 


I 

oepth 

THICK 

ALPHA 

»€TA 

9H0 

l 

S.  386*093 

5.985*903 

6.230*003 

3.497*093 

2.673*003 

2 

1.197*004 

5.985*003 

6.269*093 

3.514*993 

2.684*003 

3 

1.796*094 

5.985*003 

6.392*093 

3.S3S*0O3 

2.732*003 

4 

2.394*094 

S. 985*003 

6.522*093 

3.661*003 

2.730*093 

3 

2.993*994 

5.985*003 

6.574*093 

3.699*003 

2.753*093 

6 

3.S91*004 

5.985*003 

6.569*903 

3.637*093 

2.797*003 

7 

4 . 190*904 

5.985*003 

7.718*003 

4.332*093 

3.216*003 

9 

4.332*904 

6.417*003 

7. 7S9*003 

4.359*903 

3.227*003 

3 

5.572*004 

7.400*003 

7.328*003 

4.394*003 

3.256*003 

10 

6.429*004 

8.533*003 

7.956*003 

4.466*903 

3.303*003 

11 

7.408*004 

9.840*003 

8.101  *003 

4.5*7*003 

3.356*003 

12 

8.544*004 

1.134*004 

8.224*903 

4.616*003 

3.400*003 

13 

9.393*004 

1.308*004 

3.295*093 

4.556*993 

3.426*003 

14 

1.136*005 

1.508*904 

3.314*903 

4.567*093 

3.434*003 

IS 

1.310*005 

1.740*004 

8.397*903 

4.663*993 

3.431*003 

16 

1.511*005 

2.996*004 

8.297*993 

4.657*003 

3.427*093 

17 

1.742*005 

2.313*004 

3.298*093 

4.658*093 

3.423*003 

18 

2.008*005 

2.668*004 

8.311*003 

4.665*093 

3.432*003 

19 

2.317*095 

3.077*004 

8.329*903 

4.67S*093 

3.439*003 

29 

2.672*905 

3.S48*004 

8.347*993 

4.635*993 

3.445*003 

Figure  15.5  . 
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Figures  16.1  through  16.5  are  on  the  following  pages. 
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Figure  16.  Path  2:  SHAGAN-SHIO 

Distance:  2927  km 
Azimuth:  340.8* 

Instrument:  SRO-LP 

Events  Processed:  318,  312,  313 

Description:  This  is  an  excellent  station.  There  is  a 
small  branch  in  the  group  velocity  curve  at  0.065  Hertz, 
but  otherwise  it  is  very  clean.  The  group  velocities 
seem  stable  and  inversion  results  in  a  reasonable 
structure  with  a  relatively  low  Q. 
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SHAGAN-SHIO  STRUCTURE 


I 

DEPTH 

THICK 

ALPHA 

IETa 

OHO 

oon 

l 

5.434*903 

5.433*993 

S.  179*993 

2.397*903 

2.299*093 

9.654*991 

2 

1.087*994 

5.433*993 

5.175*993 

2.595*993 

2.238*303 

9.221*991 

3 

1 . 939*994 

5.433*993 

6.171*903 

3.464*903 

2.552*003 

1.42**002 

4 

2.174*994 

S.  433*993 

5.131*993 

3. 47S*903 

2.559*903 

2.su*9oa 

s 

2.717*994 

5.433*993 

5.237*003 

3.591*903 

2.576*903 

4.413*002 

s 

3.289*994 

5.433*993 

6.296*303 

3. =34*993 

2.537*903 

6.170*002 

7 

3.894*994 

5.433*993 

6.353*093 

3.555*993 

2.718*993 

4. 161*002 

a 

4.393*994 

5.394*993 

5.399*993 

3.532*903 

2.735*093 

2.906*002 

9 

5.974*994 

6.397*993 

6.425*993 

3.637*993 

2.746*093 

2.026*002 

t« 

5.369*994 

7.362*993 

7.732*903 

4.363*993 

3.239*993 

1.996*902 

11 

6.768*994 

9.389*993 

7.764*993 

4.355*093 

3.233*993 

1.627*002 

:2 

7.317*994 

1.348*994 

7.741*993 

4.345*993 

3.224*993 

1 . 435*002 

13 

9.928*994 

1.211*994 

7.737*993 

4.343*903 
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Figure  16.5. 
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Figures  17.1  through  17.5  are  on  the  following  pages. 


Figure  17.  Path  3:  SHAGAN-MAJO 


Distance:  4911  km 
Azimuth:  307* 

Instrument:  ASRO-LP 

Events  Processed:  318,  312,  313 

Description:  This  is  a  good  station.  There  is  clear 
evidence  of  multipathing,  but  the  separate  arrivals  are 
removed  by  the  phase  matched  filter.  The  group 
velocities  differ  by  0.1  km/sec  at  0.02  Hertz. 
Processing  three  seismograms  and  averaging  helps  to 
obtain  consistent  results. 
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IS 
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IS 
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7.S89*043 
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1.896*995 

2. *79*994 

7.S95*443 

4.263*093 

3.171*403 
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Figure  17.5. 
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Figures  18.1  through  18.5  are  on 


the  following  pages. 


0 


Figure  18.  Path  4:  SH AGAN-GRF 0 

Distance:  4699  km 
Azimuth:  62.7* 

Instrument:  SRQ-IP 
Events  Processed:  318,  312 

Description:  This  is  a  good,  although  somewhat  noisy, 
station.  The  group  velocity  plot  shows  many  distinct 
arrivals,  but  the  phase  matched  filter  seems  to  separate 
them  well.  The  final  phase  and  group  velocities  are 
quite  consistent  for  the  two  events  processed. 
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CROUP  VELOCITY  VS.  FREQUENCY 
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i 
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2.581+993 
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2 
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1.S414884 
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16 
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17 
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3.323+993 
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18 
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8.8924893 

4.542+893 
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19 
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8.1834993 

4.S93+893 

3.385+993 
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Figure  18.5 . 
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Figures  19.1  through  19.5  are  on  the  following  pages. 


Figure  19.  Path  5:  SHAGAN-ANTO 

Distance:  3740  km 
Azimuth:  57.2* 

Instrument:  SRO-LP 
Events  Processed:  318,  312 

Description:  This  is  a  low  Q  path,  so  the  signal  is 
small  relative  to  noise  and  the  data  is  difficult  to 
process.  There  are  numerous  multiple  arrivals  on  the 
group  velocity  plots.  No  group  velocity  curve  could  be 
found  for  Event  313.  Nevertheless,  fairly  consistent 
velocities  were  obtained  for  the  remaining  two  events 
from  0.02  to  0.07  Hertz,  and  when  inverted,  they  produce 
a  reasonable  velocity  and  Q  structure. 
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5.4894993 

S. 773*993 
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6 
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7 
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19 
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1 . 795*882 

Figure  19*5 


Figure  20.1  through  20.5  are  on  the  following  pages. 


Figure  20.  Path  6:  SHAGAN-CHTO 

Distance:  3890  km 
Azimuth:  337.1* 

Instrument:  SRQ-LP 

Events  Processed:  318,  312,  313 

Description:  This  is  a  rather  noisy  station  with 

multiple  arrivals  and  multipathing  present.  Most 
difficult  is  a  branch  in  the  group  velocity  curve  above 
0.06  Hertz.  This  causes  the  group  velocity  from  the 
phase  matched  filter  to  jump  to  the  other  branches. 
Averaging  velocities  for  three  events  helped  to  obtain  a 
consistent  result.  This  station  also  has  a  large  dip  in 
the  spectrum  near  20  seconds  period  which  is  difficult  to 
model  with  a  plane-layered  structure. 
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3.471*993 

2.656*093 
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3.751*993 

2.838*003 

2.241*002 

3.739*893 

2.863*003 

1.723*003 

3.877*093 

2.339*003 

1.476*093 

4.998*993 

3.005*803 

1.384*003 

4.153*993 

3.999*003 

1.397*002 
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Figure  20.5. 
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Figures  21.1  through  21.5  are  on  the  following  pages. 


Figure  21.  Path  7:  SHAGAN-KAAO 

Oi stance:  1884  km 
Azimuth:  22.1* 

Instrument:  ASRO-LP 
Events  Processed:  318 

Description:  This  is  a  very  difficult  station  to 
process  because  of  a  cleanly  bifurcated  group  velocity 
curve  above  0.04  Hertz.  The  split  is  quite  obvious  in 
the  narrow  band  filter,  but  the  phase  matched  filter 
simply  cannot  handle  it,  wanting  to  jump  back  and  forth 
between  the  branches.  To  process  this  seismogram,  we 
selected  the  upper  group  velocity  curve  and  performed  no 
iterations  with  the  phase  matched  filter.  The  final 
result  (synthetic  seismogram  made  using  the  structure) 
is  a  fairly  good  match  to  the  original,  but  we  cannot 
express  much  confidence  in  the  results.  The  same 
problem  will  arise  in  any  further  data  processing  at 
this  station;  so  it  should  probably  be  avoided  if 
possible. 
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III.  AUTOMATIC  DISCRIMINATION 


A  new  method  for  performing  seismic  discrimination  from 
teleseismic  body-wave  data  has  been  devised.  This  method  has  its 
roots  in  several  previous  methods  (VFM,  spectral  ratio  and  moment/ 
corner  frequency),  but  it  combines  and  unifies  them  into  one  master 
scheme  ideally  suited  to  automatic  analysis.  Furthermore,  it 
explicitly  incorporates  recent  results  in  deterministic  modeling  of 
explosion  and  earthquake  body  waves  and  makes  allowance  for 
absorption  through  path-dependent  t*  operators. 

For  want  of  accessible  data,  this  method,  which  we  call 
deterministic  discrimination,  has  not  been  exhaustively  tested. 
Preliminary  application  to  a  mixed  explosion/earthquake  data  set 
(Section  3.3)  shows  it  often  is  as  successful  as  the  purely 
statistical  methods.  In  order  to  apply  deterministic  discrimina¬ 
tion,  explicit  allowance  must  be  made  for  the  attenuation  of  the 
spectrum  due  to  scattering  and  absorption.  The  method  for  doing 

this,  when  the  attenuation  is  modeled  by  a  frequency-independent 
* 

t  operator  is  given  Section  3.2.  The  theory  is  then  used 
contrariwise  (Section  3.4)  to  derive  t*  directly  from  the 
deterministic  discriminant  function  evaluated  for  a  suite  of  NTS  and 
Kazakh  explosions. 

3.1  DETERMINISTIC  DISCRIMINATION 

Over  the  past  several  years  we  have  developed  the  method  for 
automatic  discrimination  outlined  in  Figure  22.  In  essence,  it 
consists  of  reducing  n  separate  recordings  of  a  seismic  event  (event 
oriented  waveform  data)  into  n  +  1  discriminant  numbers,  d,  one  for 
each  recording  as  well  as  a  network  average.  In  addition,  we  attach 
two  probabilities  to  each  d,  the  chance  d  could  have  come  from  an 
earthquake  source  p(d/q),  and  the  chance  d  could  have  come  from  an 
explosion  source  p(d/x). 

The  scalars,  d,  are  found  by  forming  a  linear  combination  of 
seismogram  features  as  is  indicated  by  the  "classify"  box  in 
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Figure  22.  Block  diagram  of  automatic  discrimination  procedure. 


Figure  22.  The  linear  combination  is  specified  by  a  set  of  weights, 
a,  and  a  constant,  b,  and  they  are  discovered  by  processing  training 
data,  as  is  shown  in  the  lower  portion  of  Figure  22. 

Heretofore,  the  training  data  has  consisted  of  seismogram 
features  calculated  several  years  ago  on  a  mixed  collection  of 
explosion  and  earthquake  data  assembled  by  VELA  Seismological  Center 
for  the  Discrimination  Experiment.  This  is  referred  to  as  the  "AI" 
data.  It  has  always  been  recognized  that  training  data  exclusively 
based  on  past  seismic  observation  has  several  deficiencies.  Among 
them  are: 

1.  For  most  collections  of  explosion  and  earthquake 

seismograms,  the  two  populations  will  have  different 
magnitude  distributions. 

2.  Rarely,  if  ever,  is  a  sufficient  number  of  events 
colocated  near  enough  that  differences  in  path  properties 
for  both  body  waves  and  surface  waves  can  be  neglected. 

3.  Methods  based  on  available  training  data  must  be  applied 
with  caution  when  analyzing  isolated  events  from  source 
regions  with  no  historical  records. 

As  a  way  of  overcoming  the  magnitude  bias,  path  bias  and 
record  absence  limitations,  and  also  to  provide  additional  informa¬ 
tion  about  the  source,  deterministic  discrimination  uses  synthetic 
training  data,  the  "Synthetic  Feature  Vectors"  shown  in  Figure  22. 

Synthetic  Feature  Vectors  describe  properties  of  the  seismic 
waves  as  they  leave  the  source  region,  not  as  they  impinge  on 
function  |a,b}  derived  from  analysis  of  synthetic  feature  vectors 
purports  to  be  universal.  This  is  in  contrast  to  our  previous 
analysis  of  AI  training  data  which  proceeded  on  a  station-by-station 
basis.  To  particularize  the  universal  discriminant  function  to  a 
specific  path,  one  geophysical  parameter,  tp,  is  needed. 

Statistical  averaging  still  is  important  for  deterministic 
discrimination,  however,  because  many  different  kinds  of  events  are 
included  in  the  synthetic  training  data.  These  events  cover  a  range 


of  source  sizes,  source  mechanisms,  and  source  material  properties. 
It  is  recognized  that  path  effects  can  blur  the  differences  between 
earthquake  and  explosion  sources,  and  an  important  part  of  the  work 
described  here  has  concentrated  on  this  point. 

Not  all  possible  discriminants  are  presently  testable  with 
synthetic  data,  for  to  create  such  data  requires  models  based  on  the 
physics  of  the  excitation  and  propagation  of  seismic  energy.  The 
complexity  discriminant,  for  example,  is  not  sufficiently  understood 
at  present  to  admit  theoretical  calculations  although  a  multiple 
rupture  earthquake  model  might  be  a  method  of  approach. 

The  synthetic  discriminant  we  have  studied  is  the  shape  of  the 
teleseismic  body  wave  spectrum  over  the  frequency  range  from  0.5  to 
5.0  hertz.  It  should  also  be  possible  to  apply  the  concept  to 

synthetic  surface  wave  spectra,  but  here  the  vagaries  of  path 

% 

dependent  attenuation  and  dispersion  require  more  particular 
information  about  specific  source-receiver  pairs  than  does  the  body 
wave  spectrum  which,  according  to  current  views,  is  controlled 
primarily  by  source  type,  source  size,  material  properties  at  the 
source  and  attenuation  (t*). 

Deterministic  discrimination  happens  in  two  stages.  The 
first,  or  learning  stage,  consists  of  training  the  algorithm  to 
recognize  the  difference  between  earthquakes  and  explosions  based 
upon  synthetic  body  wave  spectra.  The  second  stage  consists  of 
applying  the  taught  algorithm  to  actual  data  to  test  its  efficacy. 
The  only  communication  between  these  two  stages  consists  of  the 
spectral  weighting  factors  a.  Thus  we  take,  initally,  data  of  known 
type  and  find  the  weights  function  which  best  separates  their 
spectra.  These  weights  are  used  afterwards  to  classify  actual  data. 

3.1.1  Theoretical  Source  Models 

Stevens  and  Oay  (1982)  have  described  how  the  body  wave  source 
spectrums  used  in  this  analysis  were  calculated.  For  spectral 
discrimination,  of  course,  models  for  both  explosion  sources  and 
earthquake  sources  are  required.  Archambeau  et  al.  (1974)  developed 
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the  theoretical  basis  behind  spectral  discriminants  in  general. 
Analytic  approximations  to  explosion  source  spectra  have  been 
available  for  some  time  (e.g.,  Haskell,  1967;  Mueller  and  Murphy, 
1971;  von  Seggern  and  Blandford,  1972)  and  their  validity  tested 
against  actual  data  (Filson  and  Frasier,  1972;  Burdick,  et  al . , 
1981).  More  complex  explosion  models,  allowing  for  imperfect 
elasticity  and  nonlinear  stress-strain  relations,  must  be  calculated 
numerically  (Rimer,  et  aL ,  1979).  Dynamic  fault  models  of  the 
earthquake  source  are  also  a  more  recent  development,  and  likewise 


rely  on  elaborate  numerical  simulations  (e.g..  Day,  1982a,  b) . 

Numerical  models  are  expensive  to  calculate,  so  only  a  few 
distinct  source  types  are  contained  in  the  synthetic  training  data. 
Since  these  do  not  span  a  sufficiently  large  magnitude  range,  the 
few  fundamental  models  were  extended  using  cube-root  (or,  for 
explosions,  Mueller-Murphy)  scaling.  It  was  important  to  extend  the 
fundamental  source  models,  because  this  ensured  that  within  the 
training  set  there  were  sources  with  corner  frequencies  both  greater 
than  and  less  than  the  outer  limits  of  the  0.4  to  5.0  Hz  analysis 
band.  Figure  23  shows  the  steps  by  which  eight  fundamental 
explosion  models  and  six  fundamental  earthquake  models  were  extended 
to  140  synthetic  events  using  cube  root  scaling  laws.  The  final  set 
covered  a  magnitude  range  between  approximately  three  and  six. 

3.1.2  Spectral  Scaling 

It  is  a  property  of  seismic  sources  that  the  shape  as  well  as 
the  level  of  the  spectrum  changes  with  source  size.  Since  tele- 
seismic  observations  of  the  radiation  are  restricted  to  a  relatively 
narrow  frequency  band  (typically  0.5  hertz  to  3.0  hertz),  it  is 
important  to  make  a  magnitude  dependent  scaling  of  the  raw  spectra. 
The  scaling  serves  two  purposes:  it  reduces  the  influence  that 
magnitude  alone  might  have  as  a  "discriminant"  (the  magnitude  bias 
problem),  and  it  linearizes  the  difference  between  earthquake 
spectra  and  explosion  spectra,  thereby  increasing  the  efficacy  of 
the  linear  discriminant  algorithm. 


EARTHQUAKE  EXPLOSION 

MODELS  MODELS 


Figure  23.  A  handful  of  numerical  models  is  expanded  to  give  a  much 
larger  training  set  using  magnitude  scaling  laws. 


It  should  be  understood  that  the  spectral  scaling  rule  is  not 
meant  to  be  a  fit  to  any  particular  source  or  even  class  of 
sources.  It  is,  on  the  contrary,  meant  to  be  ambiguous  as  to  source 
type.  That  is,  the  rule  should  be  a  spectrum  which  is  as  much  like 
a  typical  earthquake  as  it  resembles  a  typical  explosion.  The 
fundamental  consideration  behind  the  scaling  rule  is  that  it  should 
linearize  (with  respect  to  frequency  and  source  size)  the  spectrum, 
without  showing  a  preference  for  source  type. 

The  concept  behind  the  scaling  used  here  arises  from  the  view 
that,  given  equal  moments,  ft,  an  explosion  source  has  a  larger 
corner  frequency  than  an  earthquake  source  as  shown  in  Figure  24. 
Since  the  moment  can  be  written  (Aki  and  Richards,  1980,  Equation 
14.35) 

ft  =  Df^3  ,  (1) 


it  can  be  seen  that  the  parameter,  D,  is  diagnostic  of  source  type. 
It  has,  in  fact,  been  attempted  to  infer  D  directly  from  recorded 
data  by  making  independent  measurements  of  ft  and  f  (Rivers,  et 
a!.,  1980).  The  result  was  not  particularly  successsful,  presumably 
because  of  the  limited  bandwidth  accessible  for  analysis.  Our 
approach  has  been  to  use  synthetic  modeling  in  order  to  derive  D 
from  theory  rather  than  observation;  or  more  specifically,  to 
estimate  that  D  which  best  separates  synthetic  explosion  spectra  and 
synthetic  earthquake  spectra  of  identical  moment. 

The  method  used  for  normalizing  the  spectra  is  based  on  the 
work  of  Stevens  and  Day  (1982).  In  that  study,  spectral  magnitudes 
were  found  for  a  variety  of  earthquake  and  explosion  sources  and 
used  to  evaluate  the  :MS  and  ^FM  discrimination  methods. 
Those  results  were  used  here  to  find  a  "reference  discrimination 
spectrum."  It  was  observed  that  a  spectral  magnitude  curve  of  the 
form  ( l  -  f"1) 
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with  log^D  =  5.6  cleanly  separated  the  earthquake  and  explosion 
spectral  magnitude  curves.  This  expression  was  used  to  scale  the 
theoretical  spectra  in  the  following  way.  For  each  separate 
spectrum,  the  magnitude  at  1  hertz  (mb(l))  was  noted.  This  number 
was  used  to  solve  Equation  (2)  for  the  scaling  parameter  £,  with  f 
fixed  at  1.  The  reference  spectrum,  m£(f)  was  then  evaluated 
for  all  frequencies  f  and  subtracted  from  the  actual  mb(f).  The 
spectra  used  here  varied  in  magnitude  (at  1.0  hertz)  from  3.0  to  6.0 
and  so  cover  the  important  magnitude  range.  The  resulting 

normalized  synthetic  spectra  are  shown  in  Figure  25.  Note  that  the 
normalized  explosion  spectra  almost  all  increase  with  frequency 
while  the  earthquake  spectra  generally  decrease.  Note  further,  that 
the  magnitude  zone  containing  all  spectra  of  a  given  type  is 

typically  0.8  magnitude  units  wide  at  a  given  frequency,  and  shifts, 
on  average,  by  0.5  magnitude  units  as  the  frequency  varies  from  0.5 
to  5.0  hertz.  Both  with  respect  to  the  width  of  the  spectral  band 
(^0.8  n^),  and  with  respect  to  the  magnitude  shift  with  frequency 

('v, 0-5  over  5  hz  range),  the  normalized  spectra  have  much  less 

variance  than  the  raw  spectra.  Equally  important,  the  normalization 
has  not  destroyed  the  differences  between  spectra  of  the  two 
different  source  types. 

3.1.3  Discriminating  Theoretical  Spectra 

The  scaled  theoretical  spectra  were  classified  by  using  the 
Fisher  linear  discriminant,  as  described  by  Farrell,  et  aj_.  (1981). 
A  flow  diagram  for  the  procedure  is  given  in  Figure  26. 

It  should  be  clear  from  the  scaled  spectra  shown  in  Figure  25 
that  highly  successful  discrimination  is  almost  trivially  simple. 
For  example,  when  f  <  1.0  hertz,  the  univariate  descriminant 


d  -  mb(f) 


(3) 


Figure  25.  Scaled  body  wave  spectra  for  a  collection  of 
numerical  earthquakes  and  explosions  covering 
the  magnitude  range  3  <  mb  <  6.  Dash  lines 
denote  explosions,  and  solid  lines  denote 
earthquakes. 
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Figure  26.  The  synthetic  spectra  are  scaled  and  then  classified  using  the  Fisher  discriminant. 
Misclassification  probabilities  are  estimated  by  a  jackknife  calculation. 


classifies  many  events  correctly  according  to  the  rule 
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{<  O-^explosion 

(4) 

>  0-*-  earthquake 

Simply  reversing  the  rule  gives  comparable  results  when  f  >  1.0 
hertz. 

Even  more  successful  is  the  bivariate  discriminant  defined  for 
a  pair  of  frequencies  ^f f 2*  <  **°*  *2  >  1>0} 

d  -  mb(f1)  -  mb(f2)  .  (5) 

or,  more  generally, 

d  -  a:  m^fj)  +  a2  n»b(f2)  (6) 

with  a  rule  similar  to  Equation  (4).  This  bivariate  method,  in 
fact,  was  just  the  one  developed  by  Savino,  et  a_L  (1980). 

The  multidimensional  generalization  of  the  linear  discriminant 
given  by  Equation  (6)  and  the  rule  given  by  Equation  (4)  is  the 
Fisher  discriminant 

d  ■  a  .  m  +  b  (7) 

where  in  vector  notation,  a  is  a  row  vector  of  weights,  m  is  a 
column  vector  of  normalized  magnitudes,  and  b  is  a  scalar  constant. 
To  within  an  arbitrary  scale  factor,  the  prescription  for  the  Fisher 
discriminant  says, 

a  -  S"1  mt  (8) 


a  •  (iflj  +  mg) 

b - 2 - 


(9) 


where  S  is  the  pooled  covariance  matrix,  m.  is  the  mean  vector  of 

1  t 

population  1,  m2  is  the  mean  vector  of  population  2,  and  am  is 
the  transpose  of  (m^  -  m2) . 
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In  practice,  the  pooled  covariance  is  often  nearly  singluar, 
causing  its  inverse  to  be  highly  unstable.  To  stablize  this 
calculation,  the  inverse  is  damped  (a  method  known  as  ridge 
regression  in  the  statistics  literature)  by  performing  a  singular 
value  decomposition  and  suppressing  eigenvalues  which  are  less  than 
a  certain  fraction  of  the  largest  eigenvalue.  This  leads  to  a 
slightly  nonoptimum  discriminant  function,  but  one  which  is  less 
sensitive  to  noise  in  the  data. 

The  measure  of  damping  is  tne  number  of  degrees  of  freedom 
(NDF).  When  NDF  is  small,  discrimination  is  performed  on  a  small 
dimension  subspace  of  the  space  spanned  by  all  the  variables.  For 
the  spectral  discriminant  being  studied  here,  it  is  typically  found 
that  NDF  ranging  from  2  to  10  (40  variables)  gives  satisfactory 
results,  indicating  that  the  essential  information  is  highly 
organized. 

If  it  were  reasonable  to  assume  that  the  selection  of 
synthetic  events  (the  140  spectra)  were  a  random  sample  from  two 
multidimensional  populations  of  the  Gaussian  type  with  equal 
covariances,  then  standard  formulas  could  be  used  to  estimate 
misclassification  probabilities.  We  have  not  followed  this  standard 
approach  for  several  reasons: 

1.  There  is  no  rational  justification  for  adopting  the 
multidimensional  normal  approximation  to  the  probability 
density  functions. 

2.  Even  if  there  were,  we  cannot  be  confident  that  the 
available  synthetic  data  constitutes  a  random  sample. 

3.  There  is  no  reason  to  assume  that  the  covariance  matrix 
for  explosion  spectra  is  the  same  as  the  covariance 
matrix  for  the  earthquake  spectra.  On  the  contrary,  the 
very  different  elements  of  the  source  physics  makes  it 
much  more  likely  that  collections  of  earthquake  spectra 
are  more  diverse  than  collections  of  explosion  spectra. 
Shumway  and  Blandford  (1970)  have  shown,  however,  that 
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such  differences  may  not  be  important  for  linear 
discriminant  analyses. 

Thus,  the  method  we  use  to  assign  misclassif ication  probabili¬ 
ties  is  the  jackknife  method,  a  linear  approximation  to  the 
bootstrap  (Efron,  1979,  and  references  therein).  This  technique, 
more  informatively  referred  to  as  the  leave-one-out  method,  is 
simple  to  implement.  For  a  sample  of  size  n,  where  n  =  ni  +  n2 
with  n^  samples  from  population  one  and  ng  samples  from 
population  two.  Equations  (8)  and  (9)  are  evaluated  for  all  possible 
subsets  of  the  available  training  data  containing  exactly  n  -  1 
samples.  The  set  |a,b|  is  used  in  Equation  (7),  letting  the  datum 
(m)  be  just  that  observation  deleted  from  the  |a,b(  calculation. 
Thus,  one  pretends  that  each  observation  in  turn  is  missing,  and 
classifies  it  by  the  discriminant  function  obtained  from  the 
remaining  data.  Discriminants  obtained  in  this  way,  we  refer  to  as 
d*. 

To  quantify  the  distribution  of  d*  obtained,  a  simple  normal 
law  is  fit  to  the  collection  of  each  type.  This  yields  means  and 
variances  for  the  two  populations. 

The  foregoing  procedure  was  applied  to  the  synthetic  body  wave 
spectra  yielding  the  results  plotted  in  Figure  27.  This  figure 
shows  the  jackknifed  values  for  d*  obtained  for  the  earthquakes 
(crosses)  and  explosions  (circles),  where  the  symbols  have  been 
displaced  vertically  for  clarity.  Also  shown  are  the  bell-shaped 
probability  density  functions  fit  to  each  cluster. 

The  success  of  this  procedure  is  apparent.  If  the  rule 

!<  0  explosion 

(10) 

>  0  earthquake 

is  adopted,  then  only  three  events  (one  earthquake  and  two  explo¬ 
sions)  are  misclassified.  But  these  three  events  are  among  either 
the  very  largest  or  the  very  smallest  of  the  synthetic  collection. 


90 


Figure  27.  The  Fisher  discriminant  accurately  separates  the  explosion 
spectrums  from  the  earthquake  spectrums  when  quadratic  scal¬ 
ing  (see  Figure  25)  is  used.  The  covariance  matrix  was 
highly  damped  for  this  calculation  (NDF  a  2.0).  These  jack¬ 
knife  results  show  the  explosion  discriminants  clustered 
near  y^  *  -  0.5,  and  the  earthquake  discriminants  clustered 
near  yq  ■  +  0.5.  Compare  this  figure  with  Figure  28. 


so  that  their  corner  frequencies  fall  outside  the  spectral  band  used 
for  discrimination. 

Equation  (10)  is  not  the  only  viable  rule.  It  might  be 
thought  more  hazardous  to  mistake  an  explosion  for  an  earthquake 
than  vice  versa.  In  this  case,  these  data  might  suggest  a  critical 
value  of  d*  near  +  0.25.  Alternatively,  a  rule  we  use  subsequently 
when  discussing  real  data  is  to  choose  the  critical  value  such  that 
equal  fractions  of  the  two  types  of  events  are  misclassif ied.  Of 
greater  practical  importance,  however,  is  to  discover  the  physical 
reasons  that  certain  events  are  misclassif ied,  for  this  could  lead 
to  improved  discrimination  methodologies.  The  fact  that  the  mis¬ 
classif  ied  events  shown  here  are  either  very  large  or  very  small  has 
been  noted.  Is  there,  perhaps,  a  more  subtle  residual  magnitude 
bias  hidden  in  the  distribution  of  d*?  If  so,  improvements  in  the 
spectral  scaling  law  could  be  made.  Attenuation  is  another  effect 
which  we  elaborate  on  later. 

With  the  140  test  events  contained  in  these  synthetic  training 
data,  there  are  141  possible  discriminant  functions  |a,bf,  any  or 
all  of  which  might  be  tested  against  real  data.  One  of  these  sets 
is  obtained  when  the  entire  population  of  events  is  used,  the  others 
result  from  the  separate  jackknife  calculations.  For  a  homogeneous 
collection  of  training  data,  such  as  these  synthetics,  the  differ¬ 
ences  In  the  numerous  |a,b}  are  expected  to  be  small,  particularly 
when  n  is  large.  Thus,  the  choice  ought  to  be  unimportant,  and  we 
elect  to  use  the  set  |a,bf  obtained  when  all  events  are  used.  These 
weight  factors,  called  deterministic  weights,  are  tabulated  in 
Appendix  A. 

It  was  noted  previously  that  the  judicious  choice  of  a 
reference  spectrum  is  crucial,  both  to  linearize  the  data  and  to 
ameliorate  possible  magnitude  bias.  The  quadratic  reference 
spectrum  defined  by  Equation  (2)  is  one  popular  model  which  has  some 
theoretical  justification;  but  it  is  by  no  means  the  only  reasonable 
choice.  It  is  certainly  the  simplest  model,  for  it  contains  just 
two  parameters,  D  and  f  .  However,  theoretical  studies  and 


experimental  measurements  have  shown  that  explosions  and  earthquakes 
in  some  environments  can  have  a  pronounced  spectral  rise  at  the 
corner  frequency.  This  resonant  behavior  is  evidence  of  under¬ 
damping,  an  effect  modeled  by  adding  a  term  in  the  denominator  which 
depends  linearly  on  f  .  A  cubic  denominator  polynomial  is  also 

w 

often  used,  particularly  for  modeling  earthquake  spectra. 


In  a  practical  sense,  the  necessity  of  adopting  more  complex 
reference  spectra  is,  at  present,  not  convincing.  This  is  because 
the  data  available  are  so  band  limited,  the  low  frequencies  being 
lost  in  background  noise  and  the  high  frequencies  being  lost  through 
attenuation.  For  every  extra  free  parameter  proposed  in  the 
reference  spectrum  model,  it  is  necessary  to  provide  a  means  of 
estimating  it  with  real  data. 

With  synthetic  data  this  practical  problem  does  not  occur,  and 
one  alternative  model  which  was  tested  but  not  pursued  was  the 
simple  cubic  model 


■J<f> 


lo9l0 


0f3f 

i  ♦  (*f)3 


(ii) 


The  jackknife  discrimination  result  obtained  when  this  scaling  law 

was  applied  to  the  synthetics  is  shown  in  Figure  28.  Since  the 

cubic  law  seems  to  produce  perfect  results,  why  reject  it?  First, 

the  better  performance  is  only  marginal  (3  events  in  140),  and  may 

be  accidental.  Second,  the  cubic  law  did  not  produce  better 

classification  results  than  the  quadratic  law  when  tested  on  real 

data.  Finally,  a  large  slice  of  the  mb(f)  plane  is  untouched  by 

theoretical  spectra  of  cubic  shape  (see  Aki  and  Richardo,  1980 

Figure  14.13),  and  noisy  data  falling  up  there  could,  perhaps,  be 

inadvertently  weighted  too  strongly.  The  most  mundane  reason  of  all 

is  perhaps  the  most  compelling:  real  data  is  so  band-limited,  and 

2  3 

the  corner  zone  so  broad,  that  the  asymptotic  limit  (f  or  f  ) 
is  irrelevant. 
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Figure  28.  When  cubic  scaling  Is  used,  the  jackknife  produces  this 
discriminant  distribution.  The  explosion  grouping  is 
broader  and  the  earthquake  grouping  is  narrower  than 
when  quadratic  scaling  was  used  (see  Figure  27). 


We  note  several  curious  features  on  Figure  28.  First,  the 
population  averages  and  yg  are  significantly  displaced  from 
the  reference  locations  +  0.5,  that  hold  when  the  algorithm  is  used 
on  the  entire  training  set.  Furthermore,  the  cubic  scaling  gives  a 
smaller  spread  to  the  earthquake  distribution,  but  a  larger  spread 
to  the  explosion  distribution  than  does  the  quadratic  scaling.  This 
does  not  conflict  with  the  contention  of  von  Seggern  and  Rivers 
(1979)  that  the  cubic  model  is  preferred  for  earthquakes,  but  the 
quadratic  for  explosions.  Finally,  it  should  be  noted  that  the  two 
explosions  seen  in  Figure  27  with  most  positive  d*  have  been 
gathered  into  the  main  cluster,  and  if  the  decision  rule  d*  ~  -  0.1 
is  adopted,  there  are  no  misclassif ied  events  when  the  cubic  scaling 
is  used. 

3.2  EFFECT  OF  ATTENUATION 
3.2.1  Constant  t*  Model 

When  deterministic  weights  are  used  to  classify  real  (scaled) 
body  wave  spectra  by  direct  application  of  Equation  (7),  the  two 
populations  remain  separated,  but  the  scalar  discriminant  is 
strongly  biased  in  the  positive  direction.  The  effect  is 
schematically  illustrated  by  the  idealized  plots  shown  in  Figure 
29.  The  direction  of  the  shift,  called  the  discriminant  bias,  is 
direct  evidence  that  real  spectrums  are  deficient  in  high  frequency 
energy  relative  to  the  separation  spectrum.  The  size  of  the  bias  is 
too  pronounced  to  be  attributable  to  using  a  quadratic  reference 
spectrum  in  preference  to  the  cubic.  Both  the  direction  and  size  of 
the  discriminant  bias  are  entirely  consistent  with  the  assumption 
that  it  is  a  consequence  of  frequency  dependent  attenuation  known  to 
happen  due  to  the  combined  effects  of  absorption  and  scattering  in 
the  earth. 

One  possible  way  to  handle  the  bias  that  occurs  when  determin¬ 
istic  weights  are  used  to  classify  real  data  would  be  simply  to 
measure  it  (assuming  data  from  both  populations  are  present  in  the 
training  set),  and  suitably  alter  either  the  scalar  constant,  b,  in 
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Discriminant,  d 


Figure  29.  When  deterministic  weights  are  applied  to  real  data,  the 
discriminants  are  shifted  towards  positive  values  because 
of  attenuation.  The  shift  shown  in  this  schematic  diagram 
is  exaggerated. 
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Equation  (7)  or,  equivalently,  the  threshold  value  of  d  used  in  the 
classification  rule.  We  instead  make  a  virtue  out  of  necessity  and 
show  how  the  discriminant  bias  is  not  only  a  direct  consequence  of 
attenuation,  but  can,  in  fact,  be  used  to  measure  seismic  attenua¬ 
tion  directly. 

There  are  two  points  of  view  that  may  be  adopted,  both  of 
which  lead  to  the  same  result.  From  the  receiving  station  point  of 
view,  the  cause  of  the  discriminant  bias  is  to  be  found  in  an 
incorrect  choice  of  reference  spectrum;  Equation  (2),  for  example, 
makes  no  allowance  for  the  path-dependent  spectral  attenuation  due 
to  scattering  and  absorption. 

Alternatively  (and,  being  source  modelers,  the  view  we 
prefer),  the  cause  of  the  discriminant  bias  arises  because  the 
extrapolation  of  the  far  field  signal  back  to  the  source  is 
incomplete.  Instead  of  allowing  just  for  geometric  spreading,  the 
Gutenberg  B(a)  term,  allowance  should  be  made  also  for  frequency 
dependent  attenuation.  Thus,  if  propagation  along  the  path  can  be 
modeled  by  an  operator  0(f),  then  the  VFM  magnitude  mb(f)  ought  to 
be  defined  by 

m^f)  -  log1(J[fA(f)/0(f)]  ♦  B(a)  .  (12) 


It  is  fashionable  to  posit  an  absorption  band  model  for 
explaining  seismic  attenuation,  and  Lundquist  and  Samowitz  (1982) 
have  indeed  estimated  the  parameters  in  such  a  model  from  analysis 
of  the  same  data  set  we  discuss  later.  We  have  adopted  a  more 
traditional  constant  t*  model  of  the  attenuation  function 
(Carpenter,  1967;  Oer,  e£  al_. ,  1982;  see  also  Lundquist  et  al. , 
1980,  for  a  general  review),  according  to  which 

0(f)  -  exp(-wft*)  .  (13) 

The  rule  given  by  Equation  (13)  cannot  hold  over  all 
frequencies,  but  it  is  a  good  local  approximation  over  restricted 
frequency  ranges.  The  range  here,  of  course,  is  the  teleseismic 
wave  band  0.5  to  5.0  hz. 


* 
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One  subtlety  must  be  recognized  before  Equation  (13)  can  be 
made  an  equality  and.  substituted  into  Equation  (12).  The  distance 
correction  factor  B(a)  implicitly  combines  the  effects  of  attenua¬ 
tion  as  well  as  spreading,  but  it  is  an  earth-averaged  attenuation, 
and  it  is  an  attenuation  at  the  frequency  of  1.0  hertz,  the 
predominant  frequency  present  in  the  time  domain  data  from  which 
B(a)  is  found.  To  make  this  implicit  attenuation  factor  explicit, 
we  write 

0(f)  =  exp(-ir(ft*  -  t*))  (14) 

★ 

where  t  is  the  attenuation  coefficient  for  a  particular  path 
*  P 

and  t  is  the  earth-averaged,  1.0  hertz  attenuation. 

Thus,  the  proper  way  to  define  the  variable  frequency 
magnitude  is 

n»b(f)  -  log10(fA(f))  +  1.36(f t*  -  t*)  +  B(a)  .  (15) 

It  is  the  middle  term  on  the  right  hanc  side  that  we  now  relate  to 
the  value  of  the  discriminant  bias  for  each  station. 

3.2.2  Effect  of  t*  on  the  Deterministic  Discriminant 

If  we  scale  the  spectrum  (recall  that  the  deterministic 
weights  apply  to  scaled  spectra,  not  raw  observed  spectra),  then  the 
equations  presented  in  of  Section  3.1.2  lead  to  a  modified  spectral 
magnitude  which  is  to  be  substituted  into  the  Fisher  discriminant 
given  by 

m()(f )  “  log10(fA(f))  +  BU)  +  l-36(ft*  -  t*) 

-  - 5'6  '  (16> 

In  this,  and  subsequent  equations,  the  number  5.6  is  simply 
loglo( 1 /D) *  as  discussed  in  Section  3.1.2. 
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A  problem  now  arises  —  how  is  the  source  parameter  £=f“* 
to  be  chosen?  Previously,  when  the  effect  of  attenuation  was 

I 

ignored,  i  could  be  estimated  by  requiring  mb(1.0  hertz)  =  0.0 

(see  Figure  25).  Now,  the  exact  value  of  &  is  coupled  to  assumed 
★  * 

values  for  t  and  t  .  That  good  estimates  of  source 
corner  frequencies  cannot  be  made  independently  of  estimations  of 
attenuation  is  well  known  (Filson  and  Frasier,  1972).  For  modest 
attenuation,  the  interaction  is  weak  (see  Equation  (17)  below);  so  a 
simple  iterative  analysis  method  is  adequate. 

Substituting  the  scaled  spectrum  into  the  Fisher  discriminant 
(Equation  (7)),  we  get. 


2 

d  -  a  .  [log1(J  (-  *  ^  A(f ) j  ♦  6(a)  -  5.6] 


discriminant  bias  term 


(17a) 


+  d 


bias 


(17b) 


D 

In  Equation  (17b)  we  have  defined  d  (R  for  real  data)  to  be 
the  value  of  the  discriminant  scalar  that  results  when  the  determin¬ 
istic  weights  are  applied  (incorrectly)  to  the  unmodified  VFM 
spectrum,  and  dbias  is  the  shift  caused  by  the  t*  term  in  Equation 
(17a). 

*  * 

If  tp  and  tg  are  known  quantities,  then  Equation 
(17b)  shows  how  the  initial  discriminant,  dR,  must  be  modified  so 
as  to  remove  the  shift  caused  by  attenuation.  With  reference  to 
Figure  29,  the  correction  specified  by  Equation  (17b)  takes  the  two 
Real  Data  Gaussian  functions  centered  near  3.5  (these  are  the 

D 

probability  density  functions  of  d  for  earthquakes  and  explo¬ 
sions)  and  shifts  them  leftwards  by  amount  -  dbias  to  get  two 
collections  of  d's  centered  around  zero. 


* 
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Alternatively,  if  synthetic  spectra  are  accurate  models  of 

D 

real  events,  then  observations  of  d  can  be  used  to  estimate 
tp  directly  from  the  bias  term.  This  is  the  notion  developed 
in  Section  3.4. 

One  objection  raised  during  the  initial  development  of  the  VFM 

discriminant  was  that  it  made  no  allowance  for  attenuation. 

Equation  (17b)  answers  this  objection  quantitatively.  Equation 

(17b)  also  furnishes  the  correct  framework  for  treating  data  sets 

containing  events  recorded  over  several  different  paths.  For 

example,  most  explosions  occur  in  a  few  compact  testing  areas. 

Earthquake  epicenters  are  diffuse  and  seldom  are  found  near  test 

# 

sites.  If  independent  estimates  of  tp  can  be  made  for  each 

path,  then  Equation  (17a)  can  be  used  to  equalize  the 

discriminants.  This  yields  corrected  discriminants,  d,  which  may  be 
averaged  together. 

3.2.3  Definition  of  the  Discriminant  Bias 

Before  applying  Equation  (17)  to  estimating  attenuation  from 
real  data,  it  is  necessary  to  be  more  precise  in  the  definition  of 
the  discriminant  bias.  Several  plausible  meanings  are  possible. 
Considering  first  the  synthetic  data,  define: 

C 

dj  »  the  point  where  the  explosion  probability  density 

function  p(d  X)  equals  the  earthquake  probability 
density  function  p(d  Q); 

(18a) 

S 

^2  *  the  point  where  the  explosion  misclassif ication 

probability 

(  P<dlx>) 
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equals  the  earthquake  misclassif ication 

probability 


A 

(18b) 

the  expected 

value. 

of  the  explosion 

population; 

(18c) 

the  expected 

value, 

UQ* 

of  the  earthquake 

population; 

(18d) 

(df  +  d|)/2  . 

(18e) 

These  quantities  will  be  used  in  the  left  side  of  Equation  (17b). 

Define  d^  (i  *  1,5)  for  the  real  data  in  a  similar 
manner.  Letting  DBi  denote  a  discriminant  bias,  then  the  possible 
definitions  are: 


• 

DB1 

hr  „s 

-  dl  -  dl 

(19a) 

(This  is 

the  definition  illustrated  in  Figure  29.) 

db2 

,R  hS 
■  d2  '  d2 

(19b) 

db3 

hr  ,s 

*  d3  -  d3 

(19c) 

(This  is 

the  definition  used  in  Section  3.3) 

<• 

db4 

hR  ,s 

-  d4  -  d4 

(19d) 

0B5 

,R  hS 
-  d5  -  d5 

(19e) 
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For  any  of  these  definitions,  it  then  follows  from  Equation 
(17a)  that 


DB  -  dR  -  dS 


-  1.36a  .  (ft*  -  t*)  . 


Solving  Equation  (20),  we  get 


t*  >  t*  / a  *  pB 

rp  xe  \TT7/  l.36(a  . 


T) 


where 


a  •  I  ■  Jai 


a  .  f  *2  aif. 


(20) 


(21) 


(22a) 

(22b) 


are  related  to  the  mean  and  the  first  moment,  respectively,  of  the 
deterministic  weight  coefficients.  The  expressions  given  in 
Equations  22  are  defined  over  the  whole  frequency  range,  or  any 
subrange.  Evaluating  the  expressions  with  the  deterministic  weights 
(Appendix  A)  for  two  frequency  ranges  (each  is  used  subsequently) 
gives 


a  .  I  -  -  0.1748 
a  .  f  .  -  1.425 

a  .  I  -  -  0.828 
a  .  f  «  -  3.92 


0.4  <  f  <  3.0 


0.4  <  f  <  5.0 


(23a) 


(23b) 


As  one  further  refinement,  suppose  the  discrimination 

constant,  log^D,  is  in  error  by  amount  s,  and  suppose  there  was  a 

first  approximation  to  the  path  attenuation  t  (since 

*  po 

evaluation  of  is  affected  by  choice  of  t  ,  this  method 

P 
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requires  iterative  application),  then  we  get 


l5 


t* 

po 


(24) 


A  feeling  for  the  potential  resolution  of  this  method  for 
estimating  t*  is  formed  from  consideration  of  Figure  27.  Reading 
the  jackknife  discriminant  scale  at  the  bottom  that  the  expected 
values  and  uq  are  separated  by  about  one  unit.  If  we  let 
this  be  a  measure  of  the  uncertainty  in  the  discriminant  bias,  sDB, 
then  a  first  variation  of  Equation  (24)  yields 


6tp 


«DB 

~  oirry  • 


(25) 


Using  the  first  moment  of  the  deterministic  weights  given  in 
Equation  (23),  it  follows  that 


0.4  <  f  <  3.0 
0.4  <  f  <  5.0 


(26) 


With  a  reasonable  selection  of  data  it  should  be  possible  to  do  much 
better  than  this,  particularly  if  data  of  both  classes  are  available. 


3.2.4  Classification  of  Isolated  Events 

Suppose  body  wave  spectra  are  available  for  an  event  (?)  which 
is  located  in  a  region  from  which  there  are  no  training  data. 
Suppose  we  are  unwilling  to  fix  for  certain  the  attenuation  terms 
t  for  paths  from  (?)  to  the  several  stations.  The  situation 
is  bleak,  but  not  hopeless,  for  it  is  still  possible  to  make 
quantitative  statements  about  the  likelihood  (?)  is  an  explosion  or 
an  earthquake.  This  is  done  by  using  synthetic  data,  and 
constructing  a  hypothesis  test  on  tp*. 

First,  using  any  additional  geophysical  information  such  as 
source  location  and  depth,  tectonic  regime,  surface  wave  radiation 
patterns,  etc.,  synthetic  spectra  for  a  collection  of  sources  of 


patterns,  etc.,  synthetic  spectra  for  a  collection  of  sources  of 
each  type  are  calculated.  These  synthetic  spectra  are  used  to  train 
the  discrimination  algorithm,  producing  coefficients  {a,b|  which  are 
particularized  for  the  source  region.  A  jackknife  procedure  is  used 
to  get  distribution  functions  for  d*  for  each  source  type.  These 
will  look  much  like  Figure  27.  The  deterministic  weights  are  used 
to  classify  each  spectrum  yielding  a  set  of  real  data  discrminants 

D 

d  (?)  .  Not  knowing  whether  (?)  is  an  earthquake  or  an 
explosion,  the  bias  term  cannot  be  removed,  but  we  can  get  a 
statistical  estimate  of  the  bias,  conditioned  on  whether  (?  =  X),  or 
(?  -  Q),  by  using  the  jackknife  results  obtained  on  the  synthetics. 
Thus,  for  each  observation,  either 

DB( ?  «  X)  -  dR  -  d*(X)  (27a) 

or 

DB( ?  -  Q)  -  dR  -  d*(Q)  (27b) 

since  the  d*  are  random  variables  and  so  are  the  DB.  Equations  (27) 

map  the  discriminant  bias  random  variable  onto  a  t  random 

P 

variable.  Either 

<28a> 

or 

?  - «  •  *S«ri4>-  n &V.  ?}  <28b> 

applies. 

Oeciding  whether  (?  «  X)  or  (?■  Q)  has  now  been  transformed 
into  deciding  whether  tp(  ?  «  X)  or  t*(?»Q)  is  most 
likely.  One  or  the  other,  but  not  both  are  possible,  and,  clearly, 
the  more  paths  that  are  available,  the  more  rigorous  the  test. 
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0  <  t*  <  1.5  is  a  constraint  consistent  with  almost  all  present 
P  * 

data.  Other  prior  information,  such  as  AVG(tp)  *  0.5, 

VAR(tp)  «  0.04  can  be  used  in  a  goodness  of  fit  test.  The  same 

geophysical  data  used  to  restrict  the  range  of  source  models  might 

★ 

be  used  to  restrict  the  limit  on  tp  or  to  tailor  its  mean  and 

variance.  One  example  of  this  would  be  the  correlation  between 

attenuation  and  wave  velocity  used  by  Marshall,  et  aU  (1979).  As 

attenuation  studies  mature,  we  can  anticipate  soon  the  partitioning 
* 

of  tp  into  a  source  term  and  a  receiver  term  (Der,  et  al . , 

1982b).  If  receiver  terms  are  known  for  all  stations  recording  (  ), 
then  Equation  (28)  can  be  corrected  for  them,  yielding  equations 
reflecting  just  the  source  attenuation  estimate  for  each  seismic 
recording. 

Thus,  by  rephrasing  the  problem  of  discrimination  based  upon 
teleseismic  body  waves  into  a  source  modeling,  and  t*  problem,  it  is 
possible  to  bring  into  consideration  a  variety  of  ancillary 
geophysical  data. 

3.3  APPLICATION  OF  DETERMINISTIC  DISCRIMINATION  TO  SELECTED  AI  DATA 

Deterministic  discrimination  has  been  applied  to  the  majority 
of  the  VFM  spectra  calculated  by  Savino,  ejb  aK  (1980)  during  the 
Area  of  Interest  (AI)  experiment.  To  apply  the  method  to  these 
data,  the  spectra  were  taken  exactly  as  archived:  no  new  seismic 
calculations  were  made.  Subsequent  study  of  the  AI  experiment  has 
raised  certain  questions  regarding  data  quality  and  analysis 
methodology,  but  these  we  have  not  accounted  for.  Errors  in  the 
data  are  being  corrected  and  the  data  are  being  placed  on-line.  The 
availability  of  on-line  data  will  ultimately  allow  automatic 
processing  of  scores  of  events  per  day. 

The  major  reason  for  discussing  these  results  here  is  that 
they  illustrate  the  range  of  explosion/earthquake  separation  that 
occurs  when  deterministic  discrimination  is  used  in  practice.  The 
results  range  from  excellent  to  poor.  In  some  cases  (such  as  for 
CHTO)  better  results  are  obtained  in  subsequent  analysis  (see 
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Section  3.4)  which  suggests  that  the  problems  in  the  older 
measurements  are  indeed  severe.  In  other  cases  (such  as  for  BFAK), 
the  cause  of  the  imprecision  of  deterministic  discrimination  is 
unclear. 


For  all  the  results  presented  here,  the  discriminant  function 
has  spanned  the  full  frequency  range  used  in  the  theoretical  study. 

This  is  now  known  not  to  be  optimum  for 


0.4  hertz  to  5.0  hertz. 

* 

those  paths  where  t  exceeds  approximately  0.4. 


A1 so  we 


not  performed 
calculation. 


the  iterative  improvement 


have 
★ 

required  for  best  t 

P 


3.3.1  Station  KAAO 

For  Kabul,  ten  explosion  spectrums  and  19  earthquake  spectrums 
were  used.  This  collection  corresponds  to  those  events  lying  in  the 
restricted  distance  range  15.0  <  a  <  24.0  for  which  Savino,  et  al . 
(1980)  obtained  good  separation  of  the  populations  using  the 
bivariate  (mb(high) :mb(low) )  analysis.  Equally  satisfactory 
results  are  obtained  when  the  method  of  deterministic  discrimination 
is  used:  the  two  populations  completely  separate,  i.e. ,  |d^| 
tnd  {<<§}  are  disjoint  sets.  Now,  however,  we  can  carry  tne 
analysis  a  step  farther,  turning  the  discriminant  bias  into  a  t*. 

The  VFM  (bivariate)  method  of  seismic  discrimination,  as 
originally  developed,  contained  two  empirical  factors,  the  pair  of 
fiducial  frequencies  and  the  the  equation  of  the  curve  separating 
the  two  event  populations.  The  first  subjective  factor  was  fixed  by 
scanning  scatter  plots  for  a  range  of  frequency  pairs  and  selecting 
that  pair  which  jointly  minimized  the  within-group  scatter  and 
maximized  the  between-group  scatter.  The  second  factor, 
specification  of  the  curve  separating  the  populations,  was  never 
actually  reduced  to  an  algebraic  equation,  but  decisions  were  made 
with  reference  to  an  imaginary  straight  line. 

The  deterministic  extension  of  VFM  discrimination  eliminates 
the  first  factor  (all  frequencies  are  used),  and  simplifies  the 
second  factor.  The  Fisher  discriminant,  by  projecting  the  data  onto 


106 


a  one-dimensional  subspace  (d),  rather  than  a  two-dimensional 
subspace  (m^ (high)  inflow) ),  causes  the  dimensionality  of  the 
separation  surface  to  be  reduced  from  one  to  zero.  Moreover, 
distance  perpendicular  to  the  (point)  separation  surface  has  been 
directly  related  to  a  physical  quantity,  t*. 

Using  the  second  definition  of  the  discrimination  bias  (the 
reference  point  for  each  end  of  the  line  in  Figure  29  is  equal 
misclasslfication  probability)  and  subtracting  it  from  {dR|,  gives 
the  results  shown  in  Figure  30.  Turning  the  discriminant  bias  into 
an  attenuation  (tfi  assumed  to  be  1.0)  gives  a  t*  of  0.40. 
This  represents  an  average  over  all  paths  represented  in  the  data 
set. 

Comparing  the  Kabul  explosion  data  with  the  theoretical 
explosion  data  (Figure  27),  it  is  seen  that  the  two  population  means 
are  nearly  identical  (-  0.4  vs  -  0.5,  respectively)  as  are  the 
spreads.  The  near  equality  between  explosion  means  shows  that  the 
t*  resulting  from  a  match  between  an  average  theoretical  explosion 
and  an  average  observed  explosion  (definition  DB^)  is  nearly  the 
same  as  that  obtained  using  DB2  (equal  misclassif ication 
probabilities  for  both  types  of  data)  By  Equations  (23a)  and  (23b), 
DB-j  yields  an  attenuation  of  0.42. 

It  is  somewhat  surprising  that  the  spread  of  the  synthetic 
explosion  data  matches  so  well  the  spread  of  the  observed  explosion 
data.  This  is  probably  fortuitous  for  the  synthetic  explosions 
sampled  both  a  larger  yield  range  and  a  larger  material  properties 
range  than  likely  occurs  for  the  real  data.  Further  understanding 
requires  examination  of  both  results  on  a  case-by-case  basis. 

Next  consider  the  earthquake  data.  The  assumption  that 

h 

tp  ■  0.40  shifts  the  center  of  the  observed  earthquake 
population  to  a  value  uq  -  0.85.  For  the  theoretical  earthquakes 
Mg  ■  0.50,  and  the  difference  is  negated  by  assuming 
tp  -  0.46  for  the  explosion  data.  Both  this  result  and  the 
offset  of  the  explosion  means  indicates  the  calculated  value  of  0.40 
using  definition  OBg  is  a  little  low,  but  the  difference  is  less 
than  the  error  to  be  assigned  to  the  observation. 


Figure  30.  Synthetic  weights  applied  to  KAAO  spectra  effectively 
separate  the  two  event  types.  The  discriminant  bias  is 
equivalent  to  an  attenuation  parameter  t*  ■  0.40. 
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The  most  interesting  calculation  to  be  made  is  to  test  the 
hypothesis  that  the  difference  uq  -  in  Figure  30  is  entirely 
due  to  path  differences.  It  has  been  suggested  that  this  is  a  fatal 
weakness  in  VFM  discrimination  since  real  earthquakes  and  real 
explosions  are  seldom  colocated.  To  cause  the  earthquake  distribu¬ 
tion  to  cover  the  explosion  distribution  in  Figure  30  requires  a 
shift  to  the  left  of  1.4,  corresponding  to  a  «tp  =  0.26.  Thus, 

supporting  the  hypothesis  that  the  explosion/earthquake  separation 

is  entirely  a  consequence  of  differential  attenuation  requires 
*  * 
believing  t  a  0.40  for  the  explosions  and  t  =  0.66  for 

P  *  P 

the  earthquakes.  But  the  inferred  earthquake  tp  is  a  conserva¬ 
tive  lower  bound,  for  attenuation  as  large  as  this  causes  the 
spectrum  to  fall  below  ambient  earth  noise  well  below  5.0  hertz. 

"k 

This  has  the  effect  of  biasing  tp  downwards.  This  effect  is 

accentuated  for  these  data  which  have  a  preponderance  of  small 
magnitude  earthquakes.  For  this  study  we  have  not  tested  determin¬ 
istic  discrimination  for  different  frequency  bands,  but  from  later 
calculations  we  estimate  the  bias  might  be  as  large  as  st*  =  0.15, 
pushing  the  required  earthquake  attenuation  in  the  region  of 

t*  -  0.80. 

Global  studies  of  body  wave  attenuation  do  show  extremes  in  t* 
greater  than  the  limits  (0.40  to  0.80)  required  of  the  hypothesis 
that  earthquakes  separate  from  explosions  in  these  data  because  of 
differential  attenuation.  But  analysis  of  these  events  for  other 
paths  is  required  before  a  definite  statement  can  be  made.  The 
results  in  Rivers,  et  aL  (1980)  indicate  a  much  smaller 
differential  tp. 

3.3.2  Station  CHTO 

Deterministic  discrimination  applied  to  archive  VFM  spectra 

for  the  Thailand  station  CHTO  yields  the  results  plotted  in  Figure 

31.  Using  criterion  08,  (equal  percentages  of  misclassif ied 
*  ^  * 

events)  gives  the  estimate  tp  *  0.45.  Three  out  of  11 
explosions  and  two  out  of  20  earthquakes  are  missed  if  d  =  0  is 
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taken  as  the  decision  point.  The  population  means  are  separated  by 
«d  ■  0.8  («t*  *  0.15).  The  spreads  are  significantly  larger  than 
the  spreads  found  for  either  the  synthetic  data  or  the  Kabul  data. 
For  the  explosions,  a  few  prominant  outliers  are  visible. 
Subsequent  analysis  of  a  much  larger  explosion  data  set  (see  Section 
3.4)  gives  smaller  scatter,  suggesting  some  of  these  data  are 
faulty.  This  may  also  be  the  cause  for  the  larger  scatter  in  the 
earthquake  data.  To  test  this  supposition  would  require 
case-by-case  consideration  of  these  data,  an  examination  most 
effectively  performed  interactively. 

It  should  also  be  noted  that  the  inferred  value  of  t*, 

P 

0.45,  is  considerably  greater  than  a  later  estimate  (among  results 
presented  in  Section  3.4,  is  the  value  tp  =  .25  for  CHTO), 

independent  indication  there  are  problems  in  the  data. 

It  is  concluded  that  deterministic  discrimination  for  station 
CHTO  is  less  than  perfect;  five  events  out  of  30  are  misclassified. 
Several  considerations,  however,  point  to  problems  in  the  data,  and 
re-analysis  would  probably  give  more  dramatic  results. 

3.3.3  Station  RKQN 

Red  Lake,  Ontario,  is  an  observatory  on  the  Canadian  Shield 

well  liked  by  seismologists  because  of  its  high  Q  and  impulse-like 

crustal  tr an sf erf unction.  Deterministic  discrimination  produces  the 

results  shown  in  Figure  32.  One  clear  explosion  outlier  and  one 

possible  earthquake  outlier  are  the  only  misclassified  events  out  of 

a  set  containing  six  explosions  and  38  earthquakes.  Bias  measure 

* 

DB2  produces  an  average  path  attenuation  tp  =  0.33. 

Requiring  the  five  clustered  explosion  points  to  be  centered  near 

-  0.5,  (the  mean  of  the  synthetic  explosion  data)  changes  tp  to 

0.37.  Centering  the  earthquake  population  over  the  explosion 

population  requires  a  shift  in  d  of  1.4,  equivalent  to  requiring 

* 

that  the  average  explosion  path  has  an  attenuation  constant  tp 
at  least  as  large  as  0.63.  The  bias  in  this  estimate  caused  by 
including  the  entire  0.4  to  5.0  hertz  band  is  perphaps  0.15,  as 


Figure  32.  Synthetic  weights  applied  to  RKON  spectra  produce  two 
misclassified  events,  but  the  missed  explosion  is 
clearly  anomalous.  The  parameter  t*  =  0.33. 


discussed  above.  Such  a  large  differential  attenuation  is  believed 
unlikely,  and,  as  before,  we  conclude  that  deterministic 
discrimination  is  an  effective  means  of  distinguishing  earthquakes 
from  explosions  at  this  station. 

3.3.4  Station  BFAK 

The  Alaskan  station  BFAK  yielded  the  poorest  results  (Figure 
33)  of  the  four  considered  here  but  some  of  its  neighbors  (see 
Section  3.3.5  and  Appendix  B)  do  better.  The  scatter  in  the 
explosion  data  is  particularly  pronounced  and  contrasts  sharply  with 
the  close  agreement  between  the  synthetic  data  and  the  three 
previous  examples  of  real  recordings  of  explosions. 

One  obvious  difficulty  in  analyzing  these  data  is  selecting 
the  appropriate  definition  of  the  discriminant  bias.  Choosing 
DB4,  forcing  equality  between  synthetic  earthquake  mean  and  real 
earthquake  mean,  produces  an  attenuation,  tp,  of  0.58.  This 
clearly  means  that  signals  are  falling  below  the  noise  at  the  higher 
frequencies  so  that  not  much  credence  can  be  placed  in  these  results. 

3.3.5  Other  Stations 

Deterministic  discrimination  works  well,  often  extraordinarily 
well.  This  is  remarkable  because  the  method,  as  presently  devel¬ 
oped,  has  no  adjustable  parameters.  Results  from  twelve  additional 
stations  from  the  AI  experiment  are  shown  in  Appendix  B  to  substan¬ 
tiate  this  assertion.  Not  having  settled  upon  a  particular  method 
for  forming  network  average  classification,  we  are  unable  at  present 
to  present  a  consensus  answer  for  each  event.  Furthermore,  network 
averaging  requires  regionalized  t*  information  for  best  results 
(see,  for  instance,  Der,  £t  ^1_. ,  1982b,  and  Section  3.4  of  this 
report).  More  careful  selection  of  analysis  bandwidth  and,  perhaps, 
better  spectrum  scaling  rules  are  required  as  well.  Also  of  great 
utility  would  be  an  interactive  way  of  looking  at  the  data  so  as  to 
focus  quickly  on  the  problem  events  which  should  receive  special 
scrutiny. 
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Figure  33.  Synthetic  weights  applied  to  BFAK  spectra  give  one 
mlsclasslfled  explosion  and  numerous  misclassifled 
earthquakes.  The  large  value  of  the  attenuation 
parameter  (t*  a  0.58)  is  evidence  that  too  large  a 
frequency  band  was  utilized.  Earth  noise  at  the 
higher  frequencies  may  account  for  the  relatively 
large  spread  In  these  data. 


O 


For  each  station  represented  in  the  set  of  figures  contained 

★ 

in  Appendix  B,  a  provisional  tp  estimate  is  given  based  on  the 

0  discriminant  bias.  All  data  for  a  given  station  are  plotted  with 

the  bias  removed.  Removing  a  fixed  bias,  of  course,  does  not  affect 

the  separation  between  the  earthquake  and  explosion  populations. 
★ 

The  tp  results  have  not  been  assembled  into  a  table  because  the 
^  excessive  analysis  bandwidth  and  absence  of  event  selection 

(deletion  of  deep  events,  grouping  by  source  regions)  makes  them  of 
unverified  accuracy.  (Neither  objection  applies  to  the  quantitative 
results  given  in  Section  3.4.)  The  important  relation  to  remember 
•  is  at*  »  0.19  ad,  which  aids  interpretation  of  the  plots  shown  in 

the  Appendix. 


Three  SRO  stations,  ANMO  (t*  =  0.57),  CTAO 
(t*  »  0.53),  and  ZOBO  (t*  =  0.45)  yielded  deterministic 
discriminants  (unbiased)  as  shown  in  Figures  B-l,  B-2,  and  B-3. 
These  stations  are  all  greater  than  90  degrees  from  Eastern  Kazakh, 
the  source  of  the  majority  of  the  explosions;  so  it  is  not  surpris¬ 
ing  that  the  method  is  unsuccessful.  Best  results  are  obtained  for 
CTAO,  for  if  the  deep  focus  earthquake  at  d  *  -  0.6  is  ignored 
(Event  169),  just  two  explosions  (21  and  266)  overlap  the  earthquake 
population.  We  note,  however,  the  separation  (i»q  -  vx)  is  0.5, 
only  half  the  separation  between  the  means  of  the  theoretical  data. 


Station  HNME  (Figure  B-4,  t*  -  0.51)  shows  eight 
anomalous  earthquakes  (72,  74,  24,  23,  28,  7,  9,  32),  a  well- 
clustered  earthquake  remainder,  and  a  diffuse  spread  in  the 
explosion  discriminants.  Again,  the  means  of  the  trimmed  popula¬ 
tions  have  separation  (»q  -  yx)  w  0.5. 

A  succession  of  Alaskan  stations: 


1. 

2. 

3. 

4. 

5. 


ATTU,  t 
UCAK,  t 
NJAK,  t 
TNAK,  t 
CNAK,  t 


0.29 

0.51 

0.52 

0.55 

0.49 


are  shown  in  Figures  B-5  through  B-9.  ATTU  (B-5)  looks  dismal,  but 
the  others  are  not  as  bad  as  they  seem  at  first  glance.  Earthquakes 
29  and  156  (negative  d)  as  well  as  explosions  273  and  269  (positive 
d)  are  always  grouped  in  the  wrong  class  for  these  five  stations  as 
they  were  for  BFAK.  Ignoring  those  four  events,  CNAK  ( B— 6)  exhibits 
no  errors  although  the  diffuseness  of  the  explosion  population 
suggests  two  distinct  groupings.  NJAK  ( B— 7)  is  similar  in  appear¬ 
ance,  but  a  few  events  (in  addition  to  the  four  mentioned  above)  are 
wrongly  situated.  Again,  earthquakes  are  bunched,  explosions  more 
spread  out.  For  TNAK  ( B— S )  and  UCAK  (B-9),  the  earthquake  cluster 
is  somewhat  broader  than  it  is  for  CNAK  and  NJAK,  and  the  number  of 
missed  events  correspondingly  greater. 

Three  array  sites  KSRS  (B-10,  t*  *  0.53),  Norsar  (Figure 
8-11,  t*  -  0.42)  and  LASA  (Figure  B-12,  t*  *  0.52) 
complete  the  appendix  figures.  Norsar  performs  well  although  the 
population  spreads  are  larger  here  than  they  were  for  the  synthetic 
results.  Explosion  79  (d  =  1.0),  and  earthquakes  73  and  159 
(d  ~  -  0.7)  fail  to  cluster.  There  is  some  overlap  of  the  two 
populations  near  d  -  -  0.2,  but  the  two  means  are  clearly  separated. 

The  results  shown  in  the  appendix,  in  conjunction  with  the 
four  stations  discussed  more  extensively  above,  show  that  at  many 
sites  around  the  world  deterministic  discrimination  separates  most 
earthquake  spectra  from  most  explosion  spectra.  Population 
variances  often  exceed  the  variances  in  the  theoretical  spectra, 
possibly  due  to  the  effect  of  difference  in  attenuation  along  the 
several  paths. 

If  the  average  attenuation  over  the  paths  from  all  earthquakes 
to  all  stations  exceeds  the  average  attenuation  over  the  paths  from 
all  explosions  to  all  stations  by  about  at*  =  0.15  (ad*  =  0.75), 
then  deterministic  discrimination  without  an  attenuation  correction 
is  less  reliable.  In  view  of  the  fact  that  fully  half  the  explo¬ 
sions  were  off  the  Eastern  Kazakh  test  site,  this  possibility  seems 
unlikely.  Clearly,  however,  the  more  precise  the  knowledge  of 
attenuation,  the  more  confidence  can  be  placed  on  this  technique  of 
discrimination. 
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3.3.6  Conclusions 


This  discussion  has  highlighted  significant  advances  in  the 
theory  and  application  of  seismic  discrimination  using  spectra  of 
seismic  body  waves.  Following  on  from  Stevens  and  Day  (1982),  it 
has  been  established  that  an  algorithm  trained  to  recognize  the 
differences  between  synthetic  earthquake  spectra  and  synthetic 
explosion  spectra,  with  no  modification,  is  effective  at  recognizing 
the  differences  between  real  earthquake  spectra  and  real  explosion 
spectra.  This  has  never  been  accomplished  before,  and  the  key 
ingredient  was  the  recent  development  of  dynamic  models  of  the 
earthquake  rupture  process  (Day,  1982a;  1982b).  Further  work  on 
earthquake  modeling  is  clearly  warranted. 

The  only  connection  between  the  learning  phase  of  the 

algorithm  (distinguishing  synthetic  spectra)  and  the  application 

phase  of  the  algorithm  (distinguishing  real  spectra)  was  through  a 

set  of  spectral  weights  (the  vector  a)  and  a  rule  for  spectral 
2 

scaling  (the  f  model,  D  =  5.6).  The  weights  and  scaling  law  are 

firmly  rooted  in  the  physics  of  seismic  sources  and  are  of  universal 

applicability. 

It  has  been  shown  that  the  excess  decay  of  the  spectral 
amplitudes  caused  by  attenuation  (a  mixture  of  absorption  and 

scattering)  biases  the  deterministic  discriminant.  The  bias  has 
been  interpreted  as  an  attenuation  coefficient  t*,  and  geophysically 
plausible  results  are  obtained.  Data  have  been  presented  which 
suggest,  earlier  claims  not  withstanding,  that  differential 
attenuation  does  not  account  for  the  success  of  the  VFM  discriminant. 

Deterministic  discrimination  has  been  shown  effective  under 

the  most  restrictive  of  conditions,  single  station  analysis. 
Application  of  techniques  for  combining  many  single-station  results 
into  a  network  average  can  only  improve  the  method.  This  should 
however  be  coupled  with  a  network  estimation  of  corner  frequency. 

Since  attenuation  so  quickly  drops  signals  below  noise,  an 
automatic  way  of  truncating  the  weighted  sum  (a.mb)  should  be 


developed.  This  is  a  particular  instance  of  the  more  general 
problem  of  applying  signal/noise  weighting  to  (a.m^),  and  from  the 
standard  error  in  m. (f )  at  each  frequency  calculating  an  unbiased 

0  D 

value  of  the  expected  value  of  d  and  its  standard  error. 

The  availability  of  the  universal  weights  leads  to  a  partial 
solution  of  a  long  standing  problem  in  seismic  discrimination: 
classifying  a  single  body  wave  recording  of  an  event  so  small  that 
surface  waves  are  not  seen.  The  classification  of  such  a  spectrum 
can  be  rephrased  into  a  probablility  conditioned  on  t*,  p(?  =  X/t*) 
and  p(?*Q/t*).  Refined  global  models  of  Q  in  the  crust  and  upper 
mantle  are  producing  regionalized  estimates  of  t*  of  constantly 
improved  refinement.  These  will  ineluctably  sharpen  the  method, 
even  if  only  a  general  estimate  of  source  location  is  known.  A  goal 
of  seismic  data  analysis  should  be  to  tag  each  Flinn-Engdahl  seismic 
region  with  a  t*. 

3.4  ESTIMATING  t*  FROM  SRO  RECORDINGS  OF  EXPLOSIONS 

Deterministic  discrimination  has  been  applied  to  a  collection 

of  SRO  explosion  recordings,  part  of  a  large  data  set  being 

assembled  to  study  improved  methods  of  yield  determination. 

Discrimination  is  a  misnomer  for  this  portion  of  the  study  since 

only  one  class  of  event  was  processed.  It  is  more  properly  viewed 

as  a  clustering  test  of  the  method,  for  one  purpose  was  to  see  how 

effective  the  deterministic  weights  were  at  transforming 

multidimensional  data  (VFM  spectrums)  into  compact  sets  of  numbers 

(discriminants  d).  A  second  purpose  was  to  estimate  the 

* 

discriminant  bias  and,  from  it,  infer  attenuation  parameters  tp 
on  a  station-by-station  basis  for  each  of  the  two  source  regions, 
NTS  and  Kazakh.  A  further  purpose  was  to  test  the  automatic  feature 
selection  part  of  the  program  (see  the  upper  half  of  Figure  22)  for 
the  analyses  discussed  in  Section  3.3  which  were  based  on  archive 
feature  vectors  and  involved  no  seismogram  processing. 

The  results  are  satisfying.  Nearly  500  seismograms  were 
analyzed  in  a  few  weeks.  Deterministic  weights  are  effective  at 


grouping  the  majority  of  the  data.  Attenuation  parameters,  t*,  are 
found  to  agree  well  with  other  measurements  on  these  data,  using 
entirely  different  analysis  methods. 

Some  of  the  cautions  mentioned  earlier  must  be  remembered. 

_2 

Only  the  f  spectral  model,  in  the  form  given  by  Equation  (2)  has 
been  used  to  scale  the  data.  The  corner  frequency  and  attenuation 
parameters  interact  nonlinearly,  and  we  have  not  ameliorated  this 
coupling  by  iteration.  Since  the  corner  frequency  is  a  property  of 
the  seismic  source,  it  should  be  estimated  from  all  available 
recordings  of  each  event  jointly  rather  than  on  an  individual  basis 
for  each  record.  Experience  gained  in  this  work  has  indicated  how 
many  of  these  questions  might  be  automated  with  an  improved 
processing  method. 

3.4.1  Data  Set 

The  data  set  included  38  explosions  on  the  Kazakh  test  site 
and  25  at  NTS.  The  NTS  magnitudes  for  these  events  range  from 
5*2  <  mb  <  6.2,  and  pertinent  information  is  given  in  Tables  3  and 
4.  This  is  an  important  and  much-studied  data  set,  one  of  the  most 
homogeneous  assemblages  of  digital  explosion  seismograms  yet  put 
together.  It  is  particularly  suited  to  attenuation  studies  because, 
in  addition  to  the  presumed  source  simplicity,  the  uniform  well- 
calibrated  instrumentation  makes  correction  to  earth  motion  reliable. 

Certain  deficiencies  should  be  noted,  however.  The  SRO 
stations  are  poorly  disposed  for  recording  NTS  explosions,  but  that 
is  partly  a  consequence  of  the  Pacific  Ocean  occupying  half  the 
interesting  surface  area.  Particularly  important  stations  are  KONO 
and  GRFO  in  Europe  and  MAJO  in  Japan,  for  they  pick  up  both  NTS  and 
Kazakh.  Data  from  these  three  stations  should  be  extremely  valuable 
for  differential  source  studies.  It  is  unfortunate  that  GRFO 
appears  to  lie  over  a  more  attenuative  mantle  than  normal,  and 
additional  high  latitude  stations  would  be  useful. 

The  SRO  short- period  channel  records  only  vertical  component 
motion.  Furthermore,  the  data  is  saved  only  if  the  signal  is  strong 
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Table  3 

EVENT  INFORMATION  FOR  KAZAKH  EXPLOSIONS 


Code 

Date 

mb 

SMI 

6/11/78 

5.9 

SM2 

7/5 

5.8 

SM3 

8/29 

4. 9/5. 9 

sm 

9/15 

6.0 

SM5 

11/4 

5.6 

SMS 

11/29 

5. 3/6.0 

SM7 

6/23/79 

6.3 

SMB 

7/7 

5.8 

SM9 

8/4 

6.1 

SMIO 

8/18 

6.1 

SM11 

10/28 

6.0 

SMI  2 

12/2 

6.0 

SM13 

12/23 

6.1 

SM14 

9/14/80 

6.2 

SM15 

10/29/77 

5. 5/5. 6 

SML6 

3/26/78 

5.5 

SMI  7 

4/22 

5.3 

SMI8 

7/28 

5.7 

SMI  9 

10/31 

5.2 

SM20 

5/31/79 

5.2 

SM21 

5/22/80 

5.5 

SM22 

7/31/80 

5.3 

SM23 

12/7/76 

5.9 

SM24 

5/29/77 

5.6 

SM25 

6/29/77 

5.3 

SM26 

9/5/77 

5.9 

SM27 

10/29/77 

5.6 

SM28 

11/30/77 

5.9 

Remarks 


Double  event 


Double  event 


Same  as  SM27,  double 


Table  3  (continued) 

EVENT  INFORMATION  FOR  KAZAKH  EXPLOSIONS 


Code 

Date 

mb 

Remarks 

SM29 

2/1/79 

5.4 

SM30 

4/25/80 

5.5 

SM31 

6/12/80 

5.6 

SM32 

6/29/80 

5.7 

SM33 

10/12/80 

5.9 

SM34 

12/14/80 

5.9 

SM35 

12/27/80 

5.9 

SM36 

3/29/81 

5.6 

SM37 

4/22/81 

5.9 

SM38 

5/27/81 

5.4 

<# 


<# 


* 
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Table  4 

EVENT  INFORMATION  FOR  NTS  EXPLOSIONS 


Code 

Date 

% 

Remarks 

PM25 

4/11/78 

5.3 

FONDUTTA 

PM26 

4/11/78 

5.5 

BACK BEACH 

PM27 

8/31/78 

5.6 

PANIR 

PM28 

12/16/78 

5.5 

FARM 

PM29 

6/11/79 

5.5 

PEPATO 

PM30 

9/26/79 

5.6 

SHEEP SHE AO 

PM31 

4/26/80 

5.4 

COLWICK 

PM32 

6/12/80 

5.6 

KASH 

PM33 

7/25/80 

5.5 

TAFI 

PM34 

6/6/81 

5.5 

HARZER 

YC16 

12/28/76 

5.5 

RUDDER 

YC17 

4/5/77 

5.6 

MARSILLY 

YC18 

4/27/77 

5.4 

BULKHEAD 

YC19 

5/25/77 

5.3 

CREWLINE 

YC20 

8/19/77 

5.6 

SCATLING 

YC21 

11/9/77 

5.7 

SANOREEF 

YC22 

12/14/77 

5.7 

FARALLONES 

YC23 

2/23/78 

5.6 

REBLOC HON 

YC24 

3/23/78 

5.6 

ICEBERG 

YC25 

7/12/78 

5.5 

LOWBALL 

YC26 

9/27/78 

5.7 

RUMMY 

YC27 

11/18/78 

5.1 

QUARGEL 

YC28 

2/8/79 

5.5 

QUINELLA 

YC29 

9/6/79 

5.8 

HEARTS 

YC30 

4/16/80 

5.3 

PYRAMID 
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enough  to  trigger  an  automatic  detector.  This  exacerbates  gaps  in 
coverage  due  to  maintenance  and  malfunction.  Small  magnitude  events 
are  under-represented,  not  only  because  they  often  do  not  trigger 
the  recorder,  but  also  because  it  has  been  the  recent  tendency  to 
test  in  the  intermediate  yield  range.  Even  so,  these  data  cover 
only  about  a  third  of  all  underground  shots  since  the  SRO  array  was 
deployed  starting  in  1977. 

To  show  the  general  data  quality.  Figure  34  presents  a  montage 
of  several  Kazakh  shots  with  m^  ~  5.5.  They  have  been  arranged  in 
order  of  descending  t*  (as  derived  from  the  analysis  described 
later).  The  period  lengthening,  a  consequence  of  the  low  pass 
attenuation  filter,  is  apparent,  but  the  normalized  amplitude  scale 
obscures  the  concomitant  decrease  in  amplitude.  SHIO  and  MAIO  show 
the  triplication  between  20  degrees  and  30  degrees.  Most  records 
are  simple,  but  MAJO  and  TATO  have  pronounced  codas. 

3.4.2  Feature  Extraction 

Data  were  provided  in  Lincoln  Laboratories  waveform  data  base 
format,  a  convenient  arrangement  which  collects  the  data  by  event 
and  uses  the  flexibility  of  hierarchical  directories  for  organiza¬ 
tion.  Body  wave  arrivals  were  picked  interactively  and  automati¬ 
cally  written  to  a  marker  file  for  each  event.  We  are  grateful  to 
several  Geotech  analysts  for  doing  this.  Some  events  were  repicked, 
but  the  review  was  far  from  exhaustive.  We  question  a  few  percent 
of  the  time  picks  at  most. 

Spectra  were  calculated  using  bandpass  filtering.  This  is  the 
VFM  analysis  technique  which  has  been  extensively  described  (Savino, 
et  a!.,  1980;  Farrell,  et  ah,  1981);  it  needs  little  reiteration 
here.  Each  seismogram  is  filtered  through  a  comb  of  Gaussian 
filters,  and  the  time  domain  envelope  function  found  using  the 
Hilbert  transform  method.  Amplitudes  and  arrival  times  of  envelope 
extrema  are  measured  and  corrected  for  instrument  response.  The 
ground  motion  amplitudes  are  converted  to  magnitudes  by  the  usual 
relationship. 
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CHTO  35 


BCAO  68 


ANTO  34 


SHIO  26 


MAIO  20 


MAJO  44 


TATO  41 


NWAO  90 


KONO  39 


CTAO  93 


GRFO  41 


AN  MO  95 


Figure  34.  Typical  SRO  recordings  of  Kazakh  explosions,  arranged  with 
smallest  apparent  t*  at  the  top. 


mb(f)  =  1og10(fA(f))  +  B(a) 


(29) 


ft 


ft 


ft 


ft 


ft 
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where  f  is  frequency,  A  the  earth  displacement  amplitude  (nm),  and 
B(a)  the  Gutenberg  distance  correction.  Constant  bandwidth  filters, 
Q  ®  12f  were  used  giving  a  time  resolution  of  approximately 
2.4  seconds  and  a  frequency  resolution  of  approximatey  0.13  hertz. 

The  principal  advantage  narrow  band  filtering  has  over  the 
more  usual  Fourier  analysis  is  that  a  finer  trade-off  between  time 
resolution  and  frequency  resolution  is  possible  through  the  choice 
of  filter  bandwidth.  There  is  no  linear  filter  which  gives  a 
smaller  ataf  product.  Furthermore,  it  is  not  necessary  to  engage  in 
window  carpentry  as  is  required  by  Fourier  analysis.  Ultimately  one 
probably  wants  to  tailor  the  filter  comb  for  particular  stations  and 
particular  source  regions  (paths),  but  this  has  not  yet  been  done. 
Since  the  spectra  drop  so  steeply  over  a  stoort  frequency  interval, 
side  lobe  contamination  is  a  serious  concern.  Gaussian  filters, 
however,  decay  like  exp((f  -  f  )  );  so  the  side  lobes  from  a 
signal  falling  like  exp(-f)  are  at  least  asymptotically  negligible. 
Other  things  being  equal,  side  lobe  contamination  tends  to  elevate 
the  spectral  estimates  at  high  frequencies.  This  has  the  effect  of 
biasing  t*  downwards. 

Feature  selection  proceeded  as  an  independent  process.  Once 
calculated,  VFM  spectra  were  saved  for  subsequent  use. 

The  well  recorded  signals  permitted  a  somewhat  cavalier 
attitude  towards  noise;  it  was  ignored.  Obviously  bad  data  were,  of 
course,  not  used.  Moreover,  the  exponential  decay  of  the  earth 
attenuation  operator  typically  causes  the  spectral  amplitude  to 
plunge  from  well  above  earth  noise  to  well  below  earth  noise  in  a 
fraction  of  a  hertz.  Thus,  simply  truncating  the  spectra  at  a  fixed 
high  frequency  (we  chose  3.0  hertz  and  5.0  hertz)  achieves  most  of 
the  objectives  lying  behind  more  elaborate  testing. 
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A  weighted  integral  of  scaled  spectra  in  the  form 

d  =  a  mb(f)  +  C  (30) 

was  calculated.  Since  we  are  pretending  this  was  an  assessment  of 
deterministic  discrimination,  deterministic  weights  were  taken 
(Section  3.1.3  and  Appendix  A).  The  quadratic  source  model  was  used 
for  scaling,  and  the  constant  C  was  the  discriminant  bias  term 
(Equations  (17a)  and  (17b)).  For  stations  (ANMO,  CHTO,  CTAO,  KAAO ) , 
the  specific  values  for  t*  derived  from  the  analysis  discussed  in 
Section  3.3  were  used  to  calculate  C.  For  the  remaining  SRO 
stations,  an  average  of  the  four  available  station  data  was  taken. 

Figure  35  is  typical  of  the  results  obtained.  This  figure 
plots  the  integral,  d,  vertically  against  event  number  horizont¬ 
ally.  The  dotted  vertical  lines  separate  the  results  for  KAZAKH 
(SMI  in  the  left  bin)  from  the  results  for  Yucca  Flats  (YC)  and 
Pahute  Mesa  (PM)  PMS4  in  the  right  bin).  Thus,  the  horizontal 
ordering  is  simply  by  sequence  number.  The  horizontal  dash  line 
would  separate  explosion-like  from  earthquake-like  data  if  the 
correct  value  for  t*  were  used.  This  figure  is  thus  exactly  like 
Figure  27  except  that  it  is  flipped,  turned  sideways,  and  we  have 
used  the  extra  dimension  to  stretch  out  the  individual  data  points. 
This  makes  it  easy  to  associate  the  points  with  a  particular 
seismogram. 

Comparing  again  Figures  27  and  35,  it  can  be  seen  that  we  have 
transferred  from  the  former  the  zone  covered  by  the  numerical 
explosion  models.  The  zone  is  centered  at  d  =  -  0.5. 

That  the  average  of  the  Kazakh  data  points  for  station  MAJO  is 
more  negative  than  the  average  of  synthetic  explosions  (separation 
Ad)  is  not  evidence  that  these  spectra  are  even  richer  in  high 
frequencies  than  the  numerical  models.  It  means  that  the 
provisional  correction  for  attenuation,  the  constant  C  in  Equation 
(30),  is  incorrect.  Equations  (17)  and  (21)  are  combined  to  give 
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where  d  «  -  0.5  and  d  ,  in  this  case,  is  about  -  1.1.  This 
figure  makes  obvious  the  well-known  fact  that  seismic  waves  from  NTS 
are  more  attenuated  (larger  t*)  than  seismic  waves  from  Kazakh. 
Thus,  for  each  station  we  make  separate  calculations  for  the  two 
regions. 

Analyzing  all  the  results  in  this  fashion  (Appendix  C  shows 
plots  for  the  remaining  SRO  stations)  gives  the  numbers  displayed  in 
Table  5  (Kazakh  explosions)  and  Table  6  (NTS  explosions).  In 

D 

calculating  the  averages  d  ,  a  small  fraction  of  discordant 
results  were  ignored. 

The  culling  and  grouping  has  not  been  systematic,  and  it  would 
be  extremely  valuable  to  present  displays  such  as  those  shown  in 
Figure  36  and  Appendix  C,  using  a  different  horizontal  indexsuch  as 
event  magnitude. 

Attenuation  causes  high  frequency  energy  to  be  missing  in  the 
spectrum.  Beyond  some  high  frequency  cutoff,  the  spectra  generally 
reach  the  level  of  the  to  earth  noise.  This  noise  enrichment  biases 
the  spectral  integral  negative,  yielding  an  underestimate  of  t*.  To 
explore  this  effect.  Equation  (31)  has  been  evaluated  both  for  the 
full  0.4  to  5.0  hertz  bandwidth,  as  well  as  a  more  restricted 
frequency  range  from  0.4  to  3.0  hertz.  The  figures  shown  here,  and 
the  results  in  the  Tables,  apply  to  the  restricted  interval.  (It  is 
probably  not  true  that  deterministic  weights  applied  to  synthetic 
spectra  yield  a  mean  value  d  of  -  0.5  when  only  the  3.0  hertz 
bandwidth  is  utilized.  Thus,  these  results  may  be  systematically 
biased.)  A  plot  of  t*  (0.4  -  5.0)  against  t*  (0.4  -  3.0)  is  shown 
in  Figure  36.  The  tendency  for  the  difference  between  the  t* 
produced  by  the  two  methods  to  increase  with  t*  is  evidence  of  high 
frequency  noise  contamination,  we  believe. 


Table  5 

t*  CALCULATION  FOR  KAZAKH  EVENTS 


a!*  t* 


STA 

a 

Events 

C 

u 

0 

u 

a 

ANMO 

95 

34 

1.916 

-1.104 

0.11 

0.80 

0.06 

34 

19 

1.478 

-1.413 

0.13 

0.41 

0.07 

BCAO 

68 

15 

1.478 

-1.552 

0.12 

0.34 

0.06 

BOCO 

121 

3 

1.478 

-1.343 

0.11 

0.45 

0.05 

CHTO 

35 

21 

1.275 

-1.556 

0.12 

0.24 

0.06 

CTAO 

93 

25 

1.720 

-1.188 

0.17 

0.66 

0.09 

6RF0 

41 

21 

1.478 

-0.762 

0.19 

0.75 

0.10 

GUMO 

64 

3 

1.478 

-1.030 

0.17 

0.61 

0.09 

KAAO 

18 

31 

1.000 

-0.646 

0.12 

0.56 

0.06 

KONO 

39 

20 

1.478 

-1.023 

0.30 

0.62 

0.16 

MAIO 

20 

10 

1.478 

-1.141 

0.37 

0.56 

0.19 

MAJO 

44 

25 

1.478 

-1.045 

0.18 

0.60 

0.09 

NWAO 

90 

19 

1.478 

-1.044 

0.18 

0.61 

0.09 

SHIO 

26 

13 

1.478 

-1.398 

0.15 

0.42 

0.08 

SNZO 

123 

1 

1.478 

-1.070 

0.59 

TATO 

41 

5 

1.478 

-1.098 

0.04 

0.58 

0.02 

ZOBO 

137 

21 

1.478 

-0.797 

0.15 

0.73 

0.08 

Table  6 

t*  CALCULATION  FOR  NTS  EVENTS 


STA 

A 

Events 

C 

o 

V 

a 

ANNO 

8 

20 

1.916 

-2.015 

0.10 

0.33 

0.03 

ANTO 

98 

3 

1.478 

-0.890 

0.06 

0.68 

0.03 

BCAO 

121 

5 

1.478 

-0.380 

0.20 

0.95 

0.10 

BOCO 

50 

12 

1.478 

-0.520 

0.11 

0.88 

0.05 

CHTO 

115 

6 

1.275 

-0.135 

0.18 

0.96 

0.09 

GRFO 

81 

9 

1.478 

-0.540 

0.23 

0.87 

0.12 

GUMO 

89 

2 

1.478 

-0.830 

0.01 

0.72 

0.01 

KONO 

73 

7 

1.478 

-0.544 

0.20 

0.86 

0.10 

MAJO 

79 

17 

1.478 

-0.552 

0.14 

0.86 

0.07 

ZOBO 

70 

18 

1.478 

-0.744 

0.13 

0.76 

0.07 
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Figure  36.  t*  estimates  obtained  from  body  wave  spectra  spanning  two 
frequency  Intervals. 
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In  view  of  the  confidence  limits  attached  to  the  measurements, 
these  data  may  be  grouped  as  follows: 

t*  NTS  KAZAKH 

<  0.25  CHTO 

0.30  -  0.50  ANTO,  BCAO,  SHIO 

0.55  -  0.65  CTAO,  KAAO,  KONO, 

MAIO,  MAJO,  NWAO, 

TATO 

0.70  -  0.80  ZOBO  ANMO,  GRFO 

0.85  -  1.00  BOCO,  GRFO, 

KONO,  MAJO 

Results  for  several  stations  beyond  the  distance  range  of  P  have 

been  included  (NTS  to  BCAO,  CHTO;  KAZAKH  to  BOCO,  ZOBO);  these  are 

more  useful  for  classifying  the  spectra  than  studying  the  earth's 

absorption.  The  datum  for  NTS-ANMO  is  also  suspect  because  the 
distance  is  only  eight  degrees  and  the  spectrum  may  be  distorted. 

The  tendencies  shown  in  the  above  groupings  agree  well  with 

numerous  other  studies.  More  convincing  evidence  of  the  accuracy  of 
this  procedure  can  be  seen  in  Figure  37  which  compares  our  Kazakh 
results  with  relative  attenuation  coefficients  derived  by  Lundquist 
and  Samowitz  (1982)  from  substantially  the  same  data. 

3.4.4  Conclusions 

There  are  two  important  advantages  to  this  method  of 
calculating  t*.  The  method  integrates  the  spectrum,  a  robust  and 
stable  process  and  its  uses  a  reference  spectrum  derived  in  an 
explicit  manner  from  statistically  combined  explosion  and  earthquake 
models.  Furthermore,  it  leads  to  a  compact  result,  one  number  for 
each  seismogram,  which  facilitates  the  search  for  systematic 
patterns  related  to  magnitude  bias,  depth  of  burial,  and 
regionalization.  It  is  suited  to  automatic  processing,  an  important 
consideration  when  much  data  must  be  examined. 


/ANMO 


SRO  recordings  of  KAZAKH  explosions. 


Some  weaknesses  of  the  method  have  been  mentioned,  the  major 
one,  perhaps,  being  the  use  of  single  station  corner  frequencies. 
Alternate  scaling  laws  are  possible.  It  would  be  better  to  choose  a 
somewhat  different  set  of  weight  coefficients  for  this  application, 
a  set  tuned  to  explosions  alone,  rather  than  one  tuned  to  separating 
earthquakes  for  explosions.  Also,  it  may  be  desirable  to  form  the 
synthetic  explosion  data  into  sets  based  on  similar  geologic  source 
media  (i.e.,  sets  of  explosions  in  granite,  tuff,  alluvium  and  salt) 
and  not  mix  the  sets  as  was  done  above.  Nevertheless,  both  as  a 
method  of  discrimination  and  as  a  method  of  measuring  attenuation, 
the  notion  of  a  set  of  universal,  deterministic  spectral  weights 
appears  promising. 
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APPENDIX  A 

TABLE  OF  SPECTRAL  WEIGHTS  FOR 
DETERMINISTIC  DISCRIMINATION 


APPENDIX  B 

& 

DETERMINISTIC  DISCRIMINANT  RESULTS  FOR  THE 
PRIORITY  2  AI  STATIONS  (ANMO,  CTAO,  ZOBO,  HNME, 
ATTU,  CNAK,  NJAK,  TNAK,  UCAK,  KSRS,  NAO,  AND  LAO) 
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Figure  B-2.  Deterministic  discrimination  results  for  AI  data 
recorded  at  CTAO. 
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Figure  B-3.  Deterministic  discrimination  results  for  AI  data 
recorded  at  ZOBO. 
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Figure  B-5.  Deterministic  discrimination  results  for  AI  data 
recorded  at  ATTU. 
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Deterministic  discrimination  results  for  AI  data 
recorded  at  NJAK. 
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Figure  B-10.  Deterministic  discrimination  results  for  AI  data 
recorded  at  KSRS. 


151 


-2.0  -1.5  -1.0  -0.5  0.0  0.5  1.0  1.5  2.0 

d 


Figure  8-11.  Deterministic  discrimination  results  for  AI  data 
recorded  at  NAO. 
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Figure  B-12.  Deterministic  discrimination  results  for  AI  data 
recorded  at  LAO. 


The  spectral  integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  ANMO.  Only  the  points  denoted  by  asterisks  are  meaningful. 
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Figure  C-2.  The  spectral  integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  ANTO. 
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Figure  C-3.  The  spectral  integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  BCAO. 


The  spectral  integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  BOCO. 


Figure  C-6.  The  spectral  integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  CTAO.  Only  the  points  denoted  by  asterisks  are  meaningful. 


Figure  C-7.  The  spectral  Integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  GRFO. 


The  spectral  integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  KAAO.  Only  the  points  denoted  by  asterisks  are  meaningful. 


Figure  C-10.  The  spectral  Integral  (Equation  28)  obtained  using  deterministic  weights  for 


The  spectral  Integral  (Equation  28)  obtained  using  deterministic  wei 
station  MAIO. 
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Figure  C-12.  The  spectral  integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  NWAO. 


The  spectral  Integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  SHIO. 


Figure  C-15.  ihe  spectral  integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  TATO. 
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Figure  C-16.  The  spectral  integral  (Equation  28)  obtained  using  deterministic  weights  for 
station  ZOBO.  Only  the  points  denoted  by  asterisks  are  meaningful. 


